{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating Thermodynamics Observables with a quantum computer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE, IQPE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "\n",
    "from qiskit import BasicAer\n",
    "from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.aqua.components.optimizers import SLSQP\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock\n",
    "from qiskit.chemistry.components.variational_forms import UCCSD\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "from qiskit.chemistry.core import TransformationType, QubitMappingType\n",
    "from qiskit.aqua.algorithms import VQAlgorithm, VQE, MinimumEigensolver\n",
    "from qiskit.chemistry.transformations import FermionicTransformation\n",
    "from qiskit.chemistry.drivers import PySCFDriver\n",
    "from qiskit.chemistry.algorithms.ground_state_solvers import GroundStateEigensolver\n",
    "from qiskit.chemistry.algorithms.pes_samplers.extrapolator import *\n",
    "from qiskit.chemistry.algorithms.pes_samplers.bopes_sampler import BOPESSampler\n",
    "from qiskit.chemistry.drivers.molecule import Molecule\n",
    "import qiskit.chemistry.constants as const\n",
    "from qiskit.chemistry.algorithms.pes_samplers.potentials.energy_surface_spline import EnergySurface1DSpline\n",
    "from thermodynamics_utils.thermodynamics import constant_volume_heat_capacity\n",
    "from thermodynamics_utils.vibronic_structure_fd import VibronicStructure1DFD\n",
    "from thermodynamics_utils.partition_function import DiatomicPartitionFunction\n",
    "from thermodynamics_utils.thermodynamics import Thermodynamics\n",
    "\n",
    "import warnings\n",
    "warnings.simplefilter('ignore', np.RankWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A preliminary draft with more information related to this tutorial can be found in preprint: Stober et al, arXiv 2003.02303 (2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation of the Born Oppenheimer Potential Energy Surface (BOPES) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute thermodynamic observables we begin with single point energy calculation which calculates the wavefunction and charge density and therefore  the energy of a particular arrangement of nuclei. Here we  compute the Born-Oppenheimer potential energy surface of a hydrogen molecule, as an example, which is simply the electronic energy as a function of bond length. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ft = FermionicTransformation()\n",
    "driver = PySCFDriver()\n",
    "solver = VQE(quantum_instance=\n",
    "             QuantumInstance(backend=BasicAer.get_backend('statevector_simulator')))\n",
    "me_gss = GroundStateEigensolver(ft, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))\n",
    "mol = Molecule(geometry=[('H', [0., 0., 0.]),\n",
    "                        ('H', [0., 0., 0.2])],\n",
    "                       degrees_of_freedom=[stretch1],\n",
    "                       masses=[1.6735328E-27, 1.6735328E-27]\n",
    "                       )\n",
    "\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = PySCFDriver(molecule=mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BOPES sampler testing\n",
    "bs = BOPESSampler(gss=me_gss, \n",
    "                  bootstrap=True)\n",
    "points = np.linspace(0.45, 5, 50)\n",
    "res = bs.sample(driver,points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = []\n",
    "bs_res_full = res.raw_results\n",
    "for point in points:\n",
    "    energy = bs_res_full[point]['computed_energies'][0] + bs_res_full[point]['nuclear_repulsion_energy']\n",
    "    energies.append(energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Energy')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxwUlEQVR4nO3deXxV1bn/8c83I4QpQCDMBDWoOFEN1KG1zuKIerXqr7VqbdXb+ddrq/3Z29v21lZvX51b23KtWluHWlsVlcrgbMUKKiCjDKJACAlTSBgzPL8/9ooe4gkk5Jzsk+R5v17ndfZee+29n3PE82SvtfdaMjOcc865VMqKOwDnnHNdjycX55xzKefJxTnnXMp5cnHOOZdynlycc86lnCcX55xzKefJxXV5kn4n6T9jOO8/JF19gPuOklQrKTvVcaWapB9I2iiponnckp6X9Lm4Y3QdT/6ci+vMJK0GioF6oAFYDNwHTDGzxhhDa5PwOT5nZrPijqUtJI0ClgGjzawyyfbngT+b2V0dHZuLl1+5uK7gAjPrA4wGbgduBv4Qb0hdg6Sc/VQZBWxKllhc9+bJxXUZZlZtZlOBy4GrJR0JIOleST8Iy0WSnpS0VdJmSS9Jygrbbpa0TlKNpGWSTg/l+ZJ+Lqk8vH4uKb/pvJImS5onaZuklZImhfL3m4QkHSzpWUmbQhPS/ZIKw7Y/Ef1IPxGalL4pqUSSNf24SxomaWqIeYWkzyec/7uSHpZ0X4h9kaSylr6ncNyvSFoVYvlxwndwjaR/SvqZpE3AdyX1C8eukvSupG9LypJ0BjATGBbivrd53EnO/VlJSyRtkTRd0ugD+W/tMp8nF9flmNlrwFrg40k2/0fYNoioOe3/ASbpUOBLwIRwFXQ2sDrscytwPDAeOAaYCHwbQNJEoma4bwCFwMkJ+yUS8CNgGHA4MBL4boj3KuA9oiuw3mb2P0n2fyjEPQy4FPihpNMStl8Y6hQCU4FfJzlGoouBMuBYYDLw2YRtHwVWEX0/twG/AvoBBwGfAD4DXBua8M4BykPc1+zrhJImE33flxB9/y8BD+4nTtdJeXJxXVU5MCBJeR0wlKiPoM7MXrKo47EByAfGSco1s9VmtjLs8yng+2ZWaWZVwPeAq8K264C7zWymmTWa2TozW9r8pGa2ItTZHY7xU6If6v2SNBI4CbjZzHaZ2TzgLqIf+SYvm9k0M2sA/kSUBPflDjPbbGbvAT8HrkzYVm5mvzKzemAPcAXwLTOrMbPVwE8SPn9b3Aj8yMyWhGP/EBjvVy9dkycX11UNBzYnKf8xsAKYEZqFboHoxx/4GtHVRKWkhyQNC/sMA95NOMa7oQyiK5CV7Iek4nDMdZK2AX8Gilr5WYYBm82splkMwxPWKxKWdwA99tNfsqbZsYa1sK0IyOXDnz/x3K01GvhFaJLcSvTfRwd4LJfhPLm4LkfSBKIfrJebbwt/ff+HmR1E1JT09aa+FTN7wMw+RvQjaMAdYbfyUNZkVCiD6If44FaE9cNwzKPMrC/waaIf1vdD28e+5cAASX2axbCuFedtychmxypPWE+MZSPR1V7zz38g514D3GBmhQmvnmb2ygEcy2U4Ty6uy5DUV9L5RH0Pfzazt5LUOV/SIZIEVBM1hzVKOlTSaaGjfhewE2i6lflB4NuSBkkqAr5DdOUB0V1p10o6PXRyD5d0WJLw+gC1QLWk4UR9NIk2EPVpfIiZrQFeAX4kqYeko4ma4/6crH4rfUNS/9Dk9lXgLy2cuwF4GLhNUp/QhPX1Azz374BvSToCINwocNmBhe8ynScX1xU8IamG6C/jW4n6M65toW4pMIvoh342cKeZPUfU33I70V/qFcBg4Fthnx8Ac4EFwFvAG6Gs6eaBa4GfESWrF9j7r/wm3yPqPK8GngL+3mz7j4gS2FZJNyXZ/0qghOgK41Hgv9r5TMzjwOvAvBDPvm7d/jKwnaiT/2XgAeDutp7QzB4luhp8KDQNLiS6IcB1Qf4QpXPdjCQDSkM/k3Np4VcuzjnnUs6Ti3POuZTzZjHnnHMpF8uVi6QBkmZKWh7e+7dQ7w5JC8Pr8oTy+xUNz7FQ0t2SckP5KZKqFQ3FMU/SdzrqMznnnPtALFcukv6H6KGw28NDbP3N7OZmdc4jeqjtHKI7eZ4HTjezbZLOBf4Rqj4AvGhmv5V0CnCTmZ3flniKioqspKTkwD+Qc851Q6+//vpGMxuUbNv+RjxNl8nAKWH5j0SJ4+ZmdcYRJY16oF7SAmAS8LCZTWuqJOk1YER7gikpKWHu3LntOYRzznU7kt5taVtcHfrFZrY+LFcQDZDX3HxgkqSC8ODaqez9VDGhOewq4OmE4hMkzVc0UdMRLQUg6XpJcyXNraqqateHcc45t7e0XblImgUMSbLp1sQVM7Nw3z3NymeEYTxeAaqIHnhraFbtTqKrm5fC+htEAxLWhqazx4gemvsQM5sCTAEoKyvzuxqccy6F0pZczOyMlrZJ2iBpqJmtlzQUSDrRkJndRjTkN5IeAN5OOMZ/EQ3bfUNC/W0Jy9Mk3SmpyMw2tvsDOeeca7W4msWmAk1zi19NNBTFXiRlSxoYlo8GjgZmhPXPEc23cWXiVLaShoQxo5rm2cgCNqXxczjnnEsirg7924GHJV1HNHz3JwEUzZ53o5l9jmiY75dCrtgGfDp07kM0AN67wOyw/e9m9n2iSZT+XVI90cCDV5g/yOOccx3OH6Ik6nPxu8Wcc65tJL1uZkmn1PbhX5xzzqVcXM1izrkO0tho7GlojF710auuoZH6RqO+wahvbNzrvaHRaLDwnvBqNGg0++DVGK3b++VgRO80rZth8P4y0SaMxOUPlzWxcPxEiastNbxYC3OvtVw/Q3VAy9KhQ/py3tFDU35cTy7OZZj6hkaqd9axdWcdW3fUsW1nHVt37mHrjjpqd9WzfU8D23fXs31PPdt317NjTwM79jSwc08Du+ob2F3XyK66huhV30hDY8b+dLpWkPZfpz3OO2qoJxfnOjMzo3pnHWs272TNlh2s27KTqtrdVNUkvGp3s3n7nn0eJy87i1752RTk5dA7P+f95f4FueTnZtMjJ5seuVn0yM0mPyeLvKZX9gfvudlZ5GSLnKymd5GTnUVOlshOfCl6z5LIyRZZAilaz5aQoh+/7Cwhou0IsiREU10QglC3qbzpN1Nhe+KPaFPZB8t7U0Llln57W/pRVrp/rR3gycW5lNtT38iqjbUsq6hhaUUNKytrWbNlJ2s376Bmd/1edfNzshjcN59BvfMpKSpgwpj+FPXOZ0CvPPr1zKVvz1wKe+bSL7z69MglL8e7Sl3m8+TiXDvUNTSycF01c1dv4a111SyrqGFlVS31oSkqN1uMHtiLUQMKmFjSn5EDChjRv4CRA3oyon8BfXvk+F/Srkvy5OJcG+zYU8+b723ltXc2M2f1Zt58bys766JRiYYX9uTQIX04/fDBHDqkD4cN6cuYol5+peG6JU8uzu3H1h17mLWkkqcXrufF5RvZU9+IBIcP6cvlE0YyoWQAE0r6M7hvj7hDdS5jeHJxLonKml3MWLSBpxdWMHvVJhoajWH9evB/Jo7iE4cO4rjR/enbIzfuMJ3LWJ5cnAvMjNfe2cy9r6xm+qIKGg3GFPXi+pMPYtIRQzh6RD/vH3GulTy5uG5vV10DU+eVc88rq1myfhuFBblcf/LBXPyR4Ywt7u0JxbkD4MnFdVubt+/hDy+v4sHX1rB5+x4OLe7D7ZccxeTxw+mZlx13eM51ap5cXLezp76R+2av5hfPLKd2dz1nHF7MtSeVcMJBA/0qxbkU8eTiug0z49mlldz21BJWbdzOyWMH8Z/nHU5pcZ+4Q3Ouy/Hk4rqFtzfU8N9PLual5Rs5aFAv7rlmAqccOsivVJxLk9iSi6QBwF+AEmA18Ekz25Kk3h3AeWH1v83sL6H8XuATQHXYdo2ZzQszUf4COBfYEcrfSN8ncZmsrqGRn8x4m/99aRW98rL5zvnjuOqE0eRm+4ONzqVTnFcutwDPmNntkm4J6zcnVpB0HnAsMB7IB56X9A8z2xaqfMPMHml23HOA0vD6KPDb8O66mfKtO/nyg2/y+rtbuLxsJDefcxgDeuXFHZZz3UKcyWUycEpY/iPwPM2SCzAOeDFMb1wvaQEwCXh4P8e9L0xv/KqkQklDzWx9KoN3me25ZZV8/S/z2FPfyC+v/AgXHjMs7pCc61bibBsoTvjBrwCKk9SZD0ySVCCpCDgVGJmw/TZJCyT9TFJ+KBsOrEmoszaUuW6gvqGRO55eyrX3zKG4bw+e+PLHPLE4F4O0XrlImgUMSbLp1sQVMzNJH5rRyMxmSJoAvAJUAbOBhrD5W0RJKQ+YQnTV8/02xHY9cD3AqFGjWruby2AV1bv4yoNv8trqzVw5cST/dcER9Mj151Wci0Nak4uZndHSNkkbmpqrJA0FKls4xm3AbWGfB4C3Q3nTVc9uSfcAN4X1dex9dTMilDU/7hSipERZWZlP1dfJrais5VN3vUrNrnp+fvl4LvqIX6w6F6c4m8WmAleH5auBx5tXkJQtaWBYPho4GpgR1oeGdwEXAQsTjvsZRY4Hqr2/pWtbVlHDFVNm09Bo/O3fT/TE4lwGiLND/3bgYUnXAe8CnwSQVAbcaGafA3KBl8KzCNuAT4fOfYD7JQ0imuV0HnBjKJ9GdBvyCqJbka/tkE/jYrG4fBuf/sO/yMkSD3z+eA4Z3DvukJxzgKKbqrq3srIymzt3btxhuDZasHYrV/3hNQrysnng88czpqhX3CE5161Iet3MypJt8yf0Xaf0xntbuPoPr9G3Zy4PXX88IwcUxB2Scy6BJxfX6cxZvZlr7n6Noj75PPD54xle2DPukJxzzXhycZ3K0optXHP3axT37cEDnz+eIf18amHnMpEPsOQ6ja079nD9fa/TKz/HE4tzGc6vXFynUN/QyJcffJOK6l08dIMnFucynScX1yn8ePoyXlq+kdsvOYpjR/WPOxzn3H54s5jLeI/PW8fvX1zFVceP5oqJPlSPc52BJxeX0Rauq+abjyxg4pgBfOeCcXGH45xrJU8uLmNtqt3NDX96nQG98rjzU8f6BF/OdSLe5+IyUl1DI1984A021u7mkRtPpKh3/v53cs5lDE8uLiNNeXEVr67azE8/eQxHjegXdzjOuTbydgaXcVZW1fKLZ5Zz3lFDueTYEXGH45w7AJ5cXEZpbDRu+dsCeuZm890Lj4g7HOfcAfLk4jLK/a+9x5zVW/j2eYczqI/3szjXWXlycRmjfOtObp+2hI+XFnHpcd4c5lxn5snFZQQz49uPLaTR4IcXH0WYIM4510nFklwkDZA0U9Ly8J50PA9Jd0haGF6XJ5S/JGleeJVLeiyUnyKpOmHbdzroI7l2mjq/nGeXVvIfZ431uVmc6wLiunK5BXjGzEqBZ8L6XiSdBxwLjAc+CtwkqS+AmX3czMab2XhgNvD3hF1fatpmZt9P78dwqbB5+x6+98RijhlZyLUnjYk7HOdcCsSVXCYDfwzLfwQuSlJnHPCimdWb2XZgATApsUJINqcBj6UtUpd2//3kYrbtrOOOfzuK7CxvDnOuK4gruRSb2fqwXAEUJ6kzH5gkqUBSEXAqMLJZnYuIroC2JZSdIGm+pH9IavFeVknXS5oraW5VVdWBfxLXLi+8XcWjb67jC6cewmFD+sYdjnMuRdL2hL6kWcCQJJtuTVwxM5NkzSuZ2QxJE4BXgCqi5q+GZtWuBO5KWH8DGG1mtZLOJbqiKU0Wn5lNAaYAlJWVfej8Lv0aG40fTVvC6IEFfPHUg+MOxzmXQmlLLmZ2RkvbJG2QNNTM1ksaClS2cIzbgNvCPg8AbyccowiYCFycUH9bwvI0SXdKKjKzje3+QC7lps4vZ2lFDb+4Yjz5Odlxh+OcS6G4msWmAleH5auBx5tXkJQtaWBYPho4GpiRUOVS4Ekz25WwzxCFe1glTST6fJvS8glcu+ypb+SnM9/m8KF9ueDoYXGH45xLsbgGrrwdeFjSdcC7wCcBJJUBN5rZ54Bc4KWQK7YBnzaz+oRjXBGOk+hS4N8l1QM7gSvMzJu8MtBf5q7hvc07uOeaCWR5J75zXY78tzfqc5k7d27cYXQbO/bU84kfP0/JwAIevuEEf2DSuU5K0utmVpZsmz+h7zrcva+spqpmN9+cdJgnFue6KE8urkNV76jjd8+v5LTDBjOhZEDc4Tjn0sSTi+tQv3txJdt21XPTWYfGHYpzLo08ubgOU7ltF/f88x0mjx/GuGH+wKRzXZknF9dhfvXsCuobjK+fOTbuUJxzaebJxXWIdzdt58HX3uPyCSMZPbBX3OE459LMk4vrEL98ZgU52eIrpycdjcc518V4cnFpV1G9i8fnreOKCaMo7tsj7nCccx3Ak4tLu/tmr6bBjM/6XC3OdRueXFxa7dhTzwOvvcfZ44YwaqDPMOlcd+HJxaXV395Yx9YddVz3cb9qca478eTi0qax0bj75Xc4ZkQ/ykb3jzsc51wH8uTi0ubZpZW8s3E71338IB9DzLluxpOLS5u7Xl7FsH49OOfIZBOSOue6Mk8uLi0Wrqvm1VWbueakEnKz/Z+Zc91NbP/XS7pM0iJJjWGSsJbqTZK0TNIKSbcklI+R9K9Q/hdJeaE8P6yvCNtL0vUZtu7Yw/RFFVTV7E7XKTqtu19+h1552Vw+YVTcoTjnYhDnn5QLgUuAF1uqICkb+A1wDjAOuFLSuLD5DuBnZnYIsAW4LpRfB2wJ5T8L9dLi3U07uOFPr7Ng7dZ0naJTqqjexdT55Xxywkj69cyNOxznXAxiSy5mtsTMlu2n2kRghZmtMrM9wEPAZEW9w6cBj4R6fwQuCsuTwzph++lKU29yYUH0w7llR106Dt9pNT00ee2Jfvuxc91VpjeGDwfWJKyvDWUDga1mVt+sfK99wvbqUH8vkq6XNFfS3KqqqgMKrrAgD4iax1xkx5567v+XPzTpXHeXk86DS5oFJLtV6FYzezyd594fM5sCTAEoKyuzAzlGn/wcsgTVO/3KpcnfXl9L9c46PucPTTrXraU1uZjZGe08xDpgZML6iFC2CSiUlBOuTprKE/dZKykH6Bfqp1xWlujXM5ctfuUCgJlxzz9Xc8zIQo7zhyad69YyvVlsDlAa7gzLA64AppqZAc8Bl4Z6VwNNV0JTwzph+7Ohflr0L8hjq/e5ADBn9RZWbdzOVceP9ocmnevm4rwV+WJJa4ETgKckTQ/lwyRNg/f7TL4ETAeWAA+b2aJwiJuBr0taQdSn8odQ/gdgYCj/OvD+7cvp0K8g15vFgr/MWUPv/BzOPcofmnSuu0trs9i+mNmjwKNJysuBcxPWpwHTktRbRXQ3WfPyXcBlKQ12Hwp75rKx1pvFanbVMe2t9Vz0kWEU5MX2z8o5lyEyvVks4xUW5HmfC/DkgvXsrGvgk2Uj91/ZOdfleXJpp8KCXKq9z4WH566hdHBvxo8sjDsU51wG8OTSToU986jZXU9dQ2PcocRm+YYa3nxvK58sG+kd+c45wJNLuzU9pb+tG3fqPzx3DTlZ4uJjh++/snOuW/Dk0k7dfQiYPfWN/P2NdZx++GCKeufHHY5zLkN4cmmnpiFgqnd2z079Z5dWsmn7Hi6f4B35zrkPeHJpp8Iw6m93fZDyr3PXMLhPPieXDoo7FOdcBvHk0k7duVlsw7ZdPLeskkuPG0GOTwjmnEvgvwjt1J1HRv7bG2tpNLjMn21xzjXjyaWduuvIyGbGX+euZeKYAYwp6hV3OM65DOPJpZ2aRkbubn0uc1Zv4Z2N2/2JfOdcUp5cUqB/NxwCxgepdM7tiyeXFOhuIyPX7q5n2lvrueCYoT5IpXMuqVYlF0k/kXREuoPprAq7WbPYrMUb2FnXwCXHjog7FOdchmrtlcsSYIqkf0m6UVK/dAbV2XS3kZGfmF/OsH49OG6UzzbpnEuuVcnFzO4ys5OAzwAlwAJJD0g69UBOKukySYskNUoq20e9SZKWSVoh6ZaE8vtD+UJJd0vKDeWnSKqWNC+8vnMg8bVVdxoZeeuOPby4vIrzjxlGVpYPUumcS67VfS6SsoHDwmsjMJ9oJsiHDuC8C4FLgBf3c77fAOcA44ArJY0Lm+8PcRwF9AQ+l7DrS2Y2Pry+fwCxtVl3Ghl5+qIK6hqMC44eFncozrkM1qreWEk/Ay4AngF+aGavhU13SFrW1pOa2ZJw3H1VmwisCDNOEpLYZGBxmJ2yKbbXgFgb/xNHRh7YxQdvfGL+ekoGFnDk8L5xh+Kcy2CtvXJZABxjZjckJJYmH5pqOEWGA2sS1teGsveF5rCrgKcTik+QNF/SP/Z1E4Kk6yXNlTS3qqqqXYF2lyFgqmp288rKjVxwzDCft8U5t0+tvY90PnBosx+UauBdM6tOtoOkWUCyhyBuNbPH2xRly+4EXjSzl8L6G8BoM6uVdC7wGFCabEczmwJMASgrK7P2BNFdRkb+x8L1NBpccIw3iTnn9q21yeVO4FiiKxgBRwKLgH6S/t3MZjTfwczOaGds64DEx79HhDIAJP0XMAi4IeGc2xKWp0m6U1KRmW1sZyz71F1GRn5ifjlji3sztrhP3KE45zJca5vFyoGPmFmZmR0HfARYBZwJ/E+aYpsDlEoaIykPuAKYCiDpc8DZwJVm9n4vuqQhCpdXkiYSfb5NaYrvfd2hWax8607mrN7iHfnOuVZpbXIZa2aLmlbMbDFwWFNne1tJuljSWuAE4ClJ00P5MEnTwjnqgS8B04mes3k4IYbfAcXA7Ga3HF8KLJQ0H/glcIWZtavJqzW6w8jITy1YD8D53iTmnGuF1jaLLZb0W6DptuPLQ1k+0OY/183sUeDRJOXlwLkJ69OAaUnqJY3bzH4N/Lqt8bRXdxgZ+YkF5Rw1vJ+PgOyca5XWXrlcDawAvhZeq4BriBLLAT1I2ZV09ZGRV2/czoK11VxwzNC4Q3HOdRL7vXIJDzNOM7NTgZ8kqVKb8qg6oa48MvKTC8oBOM/7W5xzrbTfKxczawAafTyxfevKIyM/MX89ZaP7M7ywZ9yhOOc6idb2udQCb0maCWxvKjSzr6Qlqk6osGcuG2u73pXLsooalm2o4XsX+qDYzrnWa21y+Xt4uRYUFuSxoqrrtRA+uaCcLME5PimYc64NWpVczOyPknoCo8yszWOJdQeFBbls3d61msXMjCfml3PCwQMZ3KdH3OE45zqR1k4WdgEwjzCGl6TxkqamMa5OpyuOjLxw3TZWb9rhD04659qstbcif5dogMqtAGY2DzgoLRF1UokjI3cV/1i4nuwscfYR3iTmnGub1iaXuiQDVHadP9FToCsOATN9UQUfHTOA/r3y4g7FOdfJtDa5LJL0f4BsSaWSfgW8ksa4Op2uNjLyisoaVlZtZ9KRftXinGu71iaXLwNHALuBB4FtRE/qu6CrjYw8fdEGAM4a58nFOdd2rb1bbAdwa3i5JJqaxbpOcqlg/MhChvTzu8Scc23X2mmOxwI3ASWJ+5jZaekJq/NpahbrCkPArNu6kwVrq7l50mFxh+Kc66Ra+xDlX4mGub8LaEhfOJ1XVxoZecaiCgDOPqI45kicc51Va5NLvZn9Nq2RdHJdaWTkpxdWMLa4NwcN6h13KM65Tqq1HfpPSPqCpKGSBjS92nNiSZdJWiSpUVLZPupNkrRM0gpJtySU3yvpnTBZ2DxJ40O5JP0y1F8g6dj2xNkWXWFk5E21u5mzejOT/NkW51w7tPbK5erw/o2EMqN9D1IuBC4Bft9ShTDc/2+IplNeC8yRNDXMhAnwDTN7pNlu5wCl4fVR4LfhPe26wsjIs5ZsoNHgLE8uzrl2aO3dYmNSfWIzWwIQprxvyURgRdN0ypIeAiYDi/exz2TgvjC98auSCiUNNbP1qYm8ZV1hZOTpizYwon9PjhjWN+5QnHOd2D6bxSR9M2H5smbbfpiuoBIMB9YkrK8NZU1uC01fPwtTLrdmHwAkXS9prqS5VVVVKQm2sCCPrZ34IcqaXXW8vHwjZx8xZH9J3znn9ml/fS5XJCx/q9m2Sfs7uKRZkhYmeU1uc6Qf9i3gMGACMAC4uS07m9kUMyszs7JBgwalIJzOPzLy88uq2NPQ6E/lO+fabX/NYmphOdn6h5jZGW2OaG/rgJEJ6yNCGQnNXLsl3UP0HM4+90m3xJGRc7Nbe69E5pi+qIKi3nkcO6p/3KE45zq5/f0CWgvLydbTYQ5QKmmMpDyiK6mpAJKGhncBFxHdIEDY/plw19jxQHVH9LdA5x4ZeVddA88treTMcUPIzvImMedc++zvyuUYSduIrlJ6hmXCervGBZF0MfArYBDwlKR5Zna2pGHAXWZ2rpnVS/oSMB3IBu42s0XhEPdLGhRimQfcGMqnAecCK4AdwLXtibMt3h8CZmcdA3vn76d2Znll5Ua272nwByedcymxz+RiZtnpOrGZPQo8mqS8nCg5NK1PI0oYzeslHXom3CX2xdRF2npNQ8Bs7YTPujy9sII++TmceHBR3KE457qAztcxkME668jI9Q2NzFpSyWmHDyYvx/9JOOfaz39JUqizjow8Z/UWNm/f4zNOOudSxpNLCnXWkZFnLK4gPyeLT4xNzS3ZzjnnySWFOuPIyGbGzMUb+NghRfTKb+1oQM45t2+eXFKoM46MvLSihrVbdnLmOL9LzDmXOp5cUiwaAqbzJJdZizcgwWmHD447FOdcF+LJJcUKC3I71a3IM5dsYPzIQgb38emMnXOp48klxQo7UbNYRfUuFqyt9iYx51zKeXJJsc40MvLMJRsAOMuTi3MuxTy5pFhnGhl51uINlAws4GCfztg5l2KeXFIscWTkTFa7u57ZKzdx5rhin7vFOZdynlxSrLOMjPxCmLvlzHH+VL5zLvU8uaRY4sjImWzWkg30L8jl2FGFcYfinOuCPLmkWGcYGbmuoZFnl1Zy2mHF5HTCSc2cc5nPf1lSrDOMjDxn9Waqd9b5LcjOubSJJblIukzSIkmNksr2UW+SpGWSVki6JaH8JUnzwqtc0mOh/BRJ1QnbvtMBH2cvnWFk5FmLK8nLyeLksT53i3MuPeIaqXAhcAnw+5YqSMoGfgOcCawF5kiaamaLzezjCfX+BjyesOtLZnZ+esLev0wfGdnMmLmkgo8dUkRBng9U6ZxLj1iuXMxsiZkt20+1icAKM1tlZnuAh4DJiRUk9QVOAx5LS6AHINNHRl62oYY1m32gSudcemVyn8twYE3C+tpQlugi4Bkz25ZQdoKk+ZL+IemINMf4IZk+MvKsxdFT+af7QJXOuTRKW7uIpFlAsocobjWzx5OUH4grgbsS1t8ARptZraRzia5oSluI73rgeoBRo0alKJxIJo+MPHOxD1TpnEu/tCUXMzujnYdYB4xMWB8RygCQVETUdHZxwjm3JSxPk3SnpCIz25gkvinAFICysjJrZ6x7ydSRkTds28X8tdV84+xD4w7FOdfFZXKz2BygVNIYSXnAFcDUhO2XAk+a2a6mAklDFMYykTSR6PNt6sCYgcwdGXmWD1TpnOsgcd2KfLGktcAJwFOSpofyYZKmAZhZPfAlYDqwBHjYzBYlHOYK4MFmh74UWChpPvBL4AozS+lVSWtk6sjIMxdvYPTAAg4Z7ANVOufSK5Z7Uc3sUeDRJOXlwLkJ69OAaS0c45QkZb8Gfp2yQA9Q1CyWWVcuNbvqeGXFJq4+cbQPVOmcS7tMbhbrtAp75lGzq576DBoZ+fkwUOXZR/hAlc659PPkkgZNT+ln0rMuMxZvoKh3Hh8Z1T/uUJxz3YAnlzTItJGRd9c38NzSSs44vJjsLG8Sc86lnyeXNMi0kZFnr9xE7e56zjrC7xJzznUMTy5pkGkjI89YvIFeedmceLAPVOmc6xieXNIgk0ZGbmw0Zi7ewCmHDqZHbnbc4TjnuglPLmnwfrNYBvS5zFu7laqa3d4k5pzrUJ5c0qBPfg652aKyZtf+K6fZ9EUV5GaLUw/zgSqdcx3Hk0saZGWJMUW9WFm5PdY4zIwZizZw/EED6dsjN9ZYnHPdiyeXNCkd3IfllTWxxrCyqpZ3Nm7nLH9w0jnXwTy5pElpcW/e27yDXXUNscUwfZEPVOmci4cnlzQpHdwHM1hRWRtbDDMWVTB+ZCHFfX3uFudcx/LkkiZji6ORh+NKLuurdzJ/bbXfJeaci4UnlzQZPbAXOVni7Q3x9Ls0TWfsA1U65+LgySVN8nKyGFPUi+UxXblMX7SBgwf14uBBPneLc67jeXJJo9Li3iyP4cqlekcdr67a5HeJOediE1tykXSZpEWSGiWV7aPe3ZIqJS1sVj5A0kxJy8N7/1AuSb+UtELSAknHpvuztKR0cJ9Y7hh7blkl9Y3md4k552IT55XLQuAS4MX91LsXmJSk/BbgGTMrBZ4J6wDnAKXhdT3w21QEeyBKi3vTaNHzJh1p+qIKivvmc8yIwg49r3PONYktuZjZEjNb1op6LwKbk2yaDPwxLP8RuCih/D6LvAoUShqagpDbbGxxHwCWb+i45LJjTz3PL6vizHHFZPncLc65mHTmPpdiM1sfliuApjag4cCahHprQ9leJF0vaa6kuVVVVWkJsCTcMdaRT+rPWlLJzroGzj96WIed0znnmktrcpE0S9LCJK/JqTyPmRlgbdxnipmVmVnZoEGDUhnO+/Jysigp6sXbHXjlMnVeOUP69mBiyYAOO6dzzjWXk86Dm9kZaTz8BklDzWx9aPaqDOXrgJEJ9UaEsliUDu7N0oqOuXKp3lHHC29XcvUJJd4k5pyLVWduFpsKXB2WrwYeTyj/TLhr7HigOqH5rMOVFvfh3U3bO+SOsacXraeuwZg8/kOtgM4516HivBX5YklrgROApyRND+XDJE1LqPcgMBs4VNJaSdeFTbcDZ0paDpwR1gGmAauAFcD/Al/okA/UgtLB0R1jq6rSP/z+1PnljCnqxZHD+6b9XM45ty9pbRbbFzN7FHg0SXk5cG7C+pUt7L8JOD1JuQFfTF2k7fP+HWOVNYwblr4f/cptu3hl5Sa+fFopkjeJOefi1ZmbxTqFkqICsrOU9tuRn1ywHjO48Bi/S8w5Fz9PLmmWn5NNycCCtA9gOXV+OeOG9uWQwT6WmHMufp5cOkDp4D5pHXr/vU07mLdmKxeO96sW51xm8OTSAcYW92Z1Gu8Ye2JBOQAXeJOYcy5DeHLpAIcU96HR4J2N6blj7PF565hQ0p/hhT3TcnznnGsrTy4doGlWynT0uyyt2MbbG2q9I985l1E8uXSAMUW9yM5SWvpdps4rJztLnHtULGNzOudcUp5cOkB+Tjaj03DHmJnxxIJyTjqkiIG981N6bOecaw9PLh2kdHDvlE95/OaarazZvNObxJxzGceTSwcZW9yHdzftYHd96u4YmzqvnLycLM4+wmecdM5lFk8uHeSQwb1paLSU3TFW39DIkwvWc/phg+nTIzclx3TOuVTx5NJBmsYYS9XcLrOWVLKxdrePgOycy0ieXDrImKJeZAlWpKBT38z43QsrGTWggDPHeZOYcy7zeHLpID1ysykZmJpZKees3sK8NVv5/MkHke2TgjnnMpAnlw50yODeLK9s/5XL719YycBeeVx23IgUROWcc6kXS3KRdJmkRZIaJZXto97dkiolLWxW/mNJSyUtkPSopMJQXiJpp6R54fW7NH+UNhlb3IfV7bxj7O0NNTyztJKrTyyhR252CqNzzrnUievKZSFwCfDifurdC0xKUj4TONLMjgbeBr6VsG2lmY0PrxtTEWyqlBZHd4yt3rjjgI8x5cVV9MzN5qrjR6cwMuecS61YkouZLTGzZa2o9yKwOUn5DDOrD6uvAp2ifah0cNMdYwfWNLa+eiePz1vH5RNG0r9XXipDc865lOoKfS6fBf6RsD5G0puSXpD08ZZ2knS9pLmS5lZVVaU/SuCgQdEdYwf6pP49/1xNo8F1HxuT4siccy61ctJ1YEmzgCFJNt1qZo+n6By3AvXA/aFoPTDKzDZJOg54TNIRZrat+b5mNgWYAlBWVmapiGd/euRmM3pgL+av2drmfat31vHAv97j/KOHMnJAQeqDc865FEpbcjGzM9J1bABJ1wDnA6ebmYVz7gZ2h+XXJa0ExgJz0xlLW1xwzDB++cxyXlmxkRMPKWr1fg/86z1qd9dz/ckHpTE655xLjU7ZLCZpEvBN4EIz25FQPkhSdlg+CCgFVsUTZXJfOOVgRg8s4NuPLWz1XWO76xu4+5/v8PHSIo4Y1i/NETrnXPvFdSvyxZLWAicAT0maHsqHSZqWUO9BYDZwqKS1kq4Lm34N9AFmNrvl+GRggaR5wCPAjWb2oRsC4tQjN5vvTz6SVRu38/sXWpf3Hn1jHVU1u7nxEwenOTrnnEsNhRalbq2srMzmzu3YlrMvPvAGMxdvYMbXTqakqFeL9RobjTN++gIF+dk88aWPIfkT+c65zCDpdTNL+qxip2wW6wq+c/448rKz+M/HF9JSgjczfvXsClZt3M4NJx/sicU512l4colJcd8e3HTWWF5avpGn3lr/oe276hr4v3+Zx89mvc2FxwzzaYydc52KJ5cYXXVCCUcN78f3n1hMza6698srt+3iiimv8ti8cm46ayy/uGK8D1DpnOtUPLnEKDtL3HbxkVTV7uYnM94GYOG6aib/5p+8vaGG3336OL50Wqk3hznnOp20PefiWufoEYVcdfxo7pu9mgG98rjz+RUM7JXPIzeeyLhhfeMOzznnDohfuWSAm84+lIG98/npzLc5clg/Hv/SSZ5YnHOdml+5ZIC+PXL51ZUf4ZWVm/jiqQeTn+ND6TvnOjdPLhni+IMGcvxBA+MOwznnUsKbxZxzzqWcJxfnnHMp58nFOedcynlycc45l3KeXJxzzqWcJxfnnHMp58nFOedcynlycc45l3I+WRggqQp4N+44UqgI2Bh3EBnCv4sP+HfxAf8uPtCe72K0mQ1KtsGTSxckaW5Ls8N1N/5dfMC/iw/4d/GBdH0X3izmnHMu5Ty5OOecSzlPLl3TlLgDyCD+XXzAv4sP+HfxgbR8F97n4pxzLuX8ysU551zKeXJxzjmXcp5cuhBJd0uqlLQw7ljiJmmkpOckLZa0SNJX444pLpJ6SHpN0vzwXXwv7pjiJClb0puSnow7lrhJWi3pLUnzJM1N6bG9z6XrkHQyUAvcZ2ZHxh1PnCQNBYaa2RuS+gCvAxeZ2eKYQ+twkgT0MrNaSbnAy8BXzezVmEOLhaSvA2VAXzM7P+544iRpNVBmZil/oNSvXLoQM3sR2Bx3HJnAzNab2RthuQZYAgyPN6p4WKQ2rOaGV7f8q1LSCOA84K64Y+nqPLm4Lk9SCfAR4F8xhxKb0BQ0D6gEZppZd/0ufg58E2iMOY5MYcAMSa9Luj6VB/bk4ro0Sb2BvwFfM7NtcccTFzNrMLPxwAhgoqRu12wq6Xyg0sxejzuWDPIxMzsWOAf4YmhaTwlPLq7LCv0LfwPuN7O/xx1PJjCzrcBzwKSYQ4nDScCFoZ/hIeA0SX+ON6R4mdm68F4JPApMTNWxPbm4Lil0Yv8BWGJmP407njhJGiSpMCz3BM4ElsYaVAzM7FtmNsLMSoArgGfN7NMxhxUbSb3CzS5I6gWcBaTsTlNPLl2IpAeB2cChktZKui7umGJ0EnAV0V+n88Lr3LiDislQ4DlJC4A5RH0u3f42XEcx8LKk+cBrwFNm9nSqDu63IjvnnEs5v3JxzjmXcp5cnHPOpZwnF+eccynnycU551zKeXJxzjmXcp5cXJckqbYVdb4mqSCF57xI0rgUHu+VduxbG96HSXpkH/UKJX3hQM/jXEs8ubju7GtAm5KLpOx9bL4ISFlyMbMTU3CMcjO7dB9VCgFPLi7lPLm4Lk3SKZKel/SIpKWS7lfkK8AwoocLnwt1z5I0W9Ibkv4axiVrmvPiDklvAJdJ+rykOWF+lL9JKpB0InAh8OPwwObBksZLelXSAkmPSuofjve8pJ9JmitpiaQJkv4uabmkHyTEXpuwfHOYd2O+pNuTfM4xIfa3mh2jpGl+H0lHhHld5oWYSoHbgYND2Y8l9Zb0TPgO3pI0OeE4SyT9b5gTZkZ42h9Jh0iaFWJ7Q9LBofwb4XtaoG4+h0y3ZGb+8leXewG14f0UoJpowMYsohEMPha2rQaKwnIR8CLRvCcANwPfSaj3zYRjD0xY/gHw5bB8L3BpwrYFwCfC8veBn4fl54E7wvJXgXKip+jzgbVNx0/4DOcArwAFYX1Aks87FfhMWP5iwr4lwMKw/CvgU2E5D+iZuD2U5xDNc9L0nawAFOrVA+PDtoeBT4flfwEXh+UeRFeDZwFTwr5ZwJPAyXH/u/BXx71yWk47znUZr5nZWoAw7HwJ0YRZiY4natL6ZzQsGXlEiajJXxKWjwxXB4VAb2B68xNK6gcUmtkLoeiPwF8TqkwN728Bi8xsfdhvFTAS2JRQ9wzgHjPbAWBmyebsOQn4t7D8J+COJHVmA7eGOU3+bmbLw2fdK3Tgh2F03EaiOXCKw7Z3zGxeWH4dKAljUw03s0dDbLvC5ziLKMG8Ger3BkqJErjrBjy5uO5gd8JyA8n/3YtozK0rWzjG9oTle4lmtZwv6Rqiq6MDjamxWXyNLcTXGvscy8nMHpD0L6LJsqZJugFY1azap4BBwHFmVhdGEO7RLGaIvsee+zidgB+Z2e/bEL/rQrzPxXVnNUCfsPwqcJKkQ+D9EWPHtrBfH2C9oiH9P5XseGZWDWyR9PGw7SrgBQ7MTODapjvbJA1IUuefRCP90iym90k6CFhlZr8EHgeOZu/vAKAf0ZwndZJOBUbvKzCLZvlcK+micI78EOd04LMJ/VbDJQ1uzYd1XYMnF9edTQGelvScmVUB1wAPhtGDZwOHtbDffxL1M/yTvYeufwj4hqQ3Q6f21UQd/AuA8UT9Lm1m0Ui1U4G5oVnvpiTVvko02dNbtDyd8yeBheEYRwL3mdkmoqbAhZJ+DNwPlIXjfIbWDc1/FfCV8DlfAYaY2QzgAWB2ONYj7J3EXBfnoyI755xLOb9ycc45l3KeXJxzzqWcJxfnnHMp58nFOedcynlycc45l3KeXJxzzqWcJxfnnHMp9/8BgjnP5rMCOwcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(points,energies)\n",
    "plt.title('Dissociation profile')\n",
    "plt.xlabel('Interatomic distance')\n",
    "plt.ylabel('Energy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_surface = EnergySurface1DSpline()\n",
    "\n",
    "xdata = res.points\n",
    "ydata = res.energies\n",
    "energy_surface.fit(xdata=xdata, ydata=ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.1576826065021262, -0.9127521770514491)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAELCAYAAAD3HtBMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr+ElEQVR4nO3deXxU1f3/8deHPQQTQBBxKy60VusO7j8fDNIKCCgqqRRatFYUxWrVoqLUVmsrUQpoEXfQr2tQKSBuLKO4s4kLWpUiKooIqCAiCMnn98fcYIAkDMncuZOZ9/PxmMfMnNzc+xmX+eScc8/nmLsjIiKyo+pFHYCIiNRNSiAiIlIjSiAiIlIjSiAiIlIjSiAiIlIjDaIOIJ1atWrl7dq1C/9CS5bAmjVw8MHhX0tEJGTz5s1b6e6tt27PqQTSrl075s6dG/6FLrwQSkogHdcSEQmZmX1cWbuGsMLQpAl8/33UUYiIhEoJJAx5ebB+fdRRiIiESgkkDE2aQGkpbNwYdSQiIqFRAglDXl7iWb0QEcliSiBhaNIk8ax5EBHJYkogYWjaNPH83XfRxiEiEiIlkDDk5yeelUBEJIspgYRBCUREcoASSBiaNUs8K4GISBZTAgmDeiAikgOUQMJQnkDWro02DhGRECmBhEE9EBHJAUogYVACEZEcoAQSBk2ii+yQ4uJi4vH4Fm3xeJzu3btv037eeedx3nnnpfzYunqO4uJiIuPuOfM44ogjPC3Kytzr13cfOjQ915OsMXz4cJ85c+YWbQMHDvSBAwdu0TZz5kzv1q3bNsdW1Z7p5xg4cKC3atVq8/EzZ870Vq1a+YgRI7ZpLygo8MLCwpQfW1fPsfU/4zAAc72S79TIv9TT+UhbAnF3Lyhwv/ji9F1PMlJlCaG6L93KvkhT9uVTUODxadPc163zF6ZO9X1atvQxf/ub/7xlS3+ppMT9k0/85Yce8iNatPB7rrrKj2zRwl8bP9594UJ//Z57/Oj8fD+uWTOffccd7vPm+ezbb/cTCwv9gYsv9l8WFvrcMWPcX33V5956q3crLPSHBw/27oWFPm/0aPdZs3zeqFHetWlT756f7/NHjnSPx33+v/7lvQoKfP6IET5/xAg/paDA7+3f308tKPD5N9/sPmOGz7/5Zj+1oMDH9e/vvQsK/I2bb/Y3br7ZewdtpxUU+Bs33eQ+bZq/cdNNflpBgY/r129ze2Vt5ceeXlDg4/v189MLCvyN4mL3557zN4qLt2mvrK2qY6M4xxaPadMqf6xcWeP/jpVA0p1A2rZ1P+ec9F1PIldZshgxYoTn5+dv90s+/swzfkDLlv7auHE+95Zb/NcFBf5o797+p2bN/MNBg3zxgAF+a16ev37EEf5w48a+/IQT3Hv08FWHH+6vNWjgS9u29YX16/t3u+/uvtdevr5lS19l5t83auTrwEvr1Uv8765H7j6efrrG/20rgXiaE8h++7mfeWb6ridpU1WvorLewx477+z3DhnipxcU+MQePfzGpk19ac+e7r16+Tc//7kvrlfPv2/UKKkvgA0NGviX4F8VFrrvv7/7YYe5H3usL9p7b58C/s7++7v37et+1lnuAwf6ax06+L/AXzzmmMRw6rXX+vROnfxK8Ge6dHEfOdL9llvcx4zxSd27+zngT/Ts6T5+vPv997s/+KA/2ru3F4E/cvrp7o895v744/5gnz7eC/yBoiL3yZPdp0xxf/JJv//MM70r+H19+7o/84z7s8+6P/ecj+/Xz08EH9e/v/v06e4zZvi9v/2tdwK/57e/dY/HE72RkSO9V2Gh3/O733nPwkKfN3Kk+wsv+LyRI71nYaHf/bvfeY/CQp83apTPGzXKexQW+t0DBiTaRo92f/FFnzd6tJ9cWOh3DRjgJwftlbVtcexZZyV6Srfc4v7SSz7vllu8+1btlbVVdWwU59j8ePHFqh+rVtX4v3klkHQnkEMOce/ZM33Xk5RLOlHMmOEHtWjhC/75T/9w0CD/v8aN/aO99vLPq/qrv1Ur94MOcu/SxRf84hc+Enx6p07uY8a4P/ywL7jxRj+psNBHn3++/6JFC39hyhSfOW2at2rVyocNG1bpsFUy7TtybLrPsfV4vuZANAeScY+0JpDjjnPv3Dl915OUq+pL7cUnnvA3r7/eR+Xl+aK99/YVZlskiG+bNvVZ4PMPPtj9uut84VVXec/CQh9x0UXeduedk/pyrUtfYKk4x8CBA5OeK9qRyfxMv3kgFecYPny4h00JJN0J5Fe/cj/qqPRdT2qlut7GPi1b+kN9+vg9TZr42j333JwoNtWr53PB5xx2mPu//+0+a5bPmjhxi6SwI1+uVX2RZvoXWF35EpSaUwJJdwLp3dv9wAPTdz2pla3/8n/50Uf9yvx8X3X44b4pGIpa37Che/fu7sOH+7zRo32PnXfebu+hadOmPmLEiG2uVdWXrr5IJRMpgaQ7gfTv796uXfquJ7X2wpQpflGzZv6/du28NOhlrN1rLx+dl+d3DxiwefipqqGtqoZhlBSkrqsqgWgleljy87USPQNVtuL59fHjmX/ssZzQty+3rF1L2ZIlPH/CCbw+fjzt1q3joKlTOWf8eB6cMIGioiIeeeQRSkpKiMViAMRiMUpKSth33303t5WLxWIMGTIkbZ9PJJ0aRB1A1mrWTAkkA3Xs2JGioqJEAmjYkK/+9CeOmjuX0oYNWda5M79/7TU6XnQRY2+/ndNeeaXSRDFnzpxKE8XWbSLZTgkkLPn5sG4dlJVBPXX0MkUsFuPZa67hu5NOgo0b2WTG4nPO4fOTT6b3wIGUTJyYSAadO1NUVMSZZ565ze8rUYgk6JstLOUVeb//Pto4clRlQ1WvPvQQ7x90EIdfcgmHNm7MJcAdV1zBPnffzSsfflhlb0NEKqcEEhZtKhWp8qGqeDwO69ezpH9/DuvXj/0WLWLx2WdzQOPGFAwbxi133008HmfIkCGavxDZQZEkEDNraWbTzOzD4LlFFccNN7N3gsevK7SPN7OPzGxB8Dg0bcEnS3uCRKq8B/HP3r35cs89affgg3zTpQuv3XcfR02Zwv0TJnDddddRUlLyY6IRkR0SVQ/kSmCGu7cHZgTvt2BmJwOHA4cCRwGXm1lBhUP+7O6HBo8F4Ye8g5RAorVxI7Fp03hmzRo2rFzJ/X37suu0aby8ZImGqkRSJKpJ9FOATsHr+4DngSu2OuYAYJa7bwI2mdlbQFegJE0x1o42lUqb4uJiOnbs+OMQ1Oef803XrjR/+20ebtyYJX/8I6PGjWPPYKhqa5oYF6mZqHogbdx9WfD6C6BNJce8CXQ1s6Zm1gqIAXtW+PkNZvaWmY00s8ZVXcjMBprZXDObu2LFipR9gO1SDyRttpjveOEFNvziFzR8+23+kJfHbk8/zdXFxRqqEglBaAnEzKZXmL+o+Dil4nHBKkff+vfd/TngKeAV4GHgVaA0+PFVwP5AR6Al2/ZeKp7nTnfv4O4dWrdunZLPlhRNoqdN+TDUxF69KI3F+HTNGm464wz6TZ2qoSqREIU2hOXuXar6mZktN7O27r7MzNoCX1ZxjhuAG4LfeQj4IGgv771sMLNxwOUpDT4V1ANJH3diL79MbO1angPmXHYZfx0+fJvDNFQlklpRDWFNBgYErwcAk7Y+wMzqm9nOweuDgYOB54L3bYNnA04F3gk/5B2kBJIepaUwaBAMG8ajjRvz6tChjLr3Xg1ViaRBVJPoNwIlZnYO8DFQBGBmHYDz3f0PQEPgxUSOYA3QP5hQB3jQzFoDBiwAzk9v+EnQJHr4ysrg3HNh3DhG5+Vx8JNP8uvOnTmhS5cfy5WoxyESmkgSiLuvAk6spH0u8Ifg9XoSd2JV9vudQw0wFdQDSbkt7rYqK4Pzz4dx43hwv/04+M47k6pZJSKpo1pYYWnUCBo00CR6Cm0uhPjoo8QefxzuuouReXkcescdKm4oEgElkDCppHtKlfcsXuvRg9i6ddyal8ehTz5JrHPmd0hFspESSJiUQFIu9vnnxNat4z5gxWWXKXmIREjFFMOkBJJaL75I2dln81LDhiwZOpSxt9+uu61EIqQEEiZtKpU6H33Exh49+F9ZGT5hAtfecINWl4tETAkkTOqBpMb69dCnD6UbN7Jy/Hj+3ymJYgZaXS4SLSWQMOXn6y6sGtpiQ6jLLoN58/hg6FBe/PzzLY7Tnh0i0VECCZN6IDVWfsvuwmHD4Lbb+KSoiBNHj6Zjx45RhyYiAd2FFSYlkBqLxWJMGTmSn/zud3yyxx4cNWMGJRMmaG2HSAZRDyRMmkSvudJSjr79duo3asRxS5dy7gUXKHmIZBglkDCpB1Jzo0fDyy9zecOGnD1sGGPHjtXdViIZRgkkTPn58P33iYqxkrz//pfSK6/k6UaNOGPSJO1dLpKhlEDCVF5Qcd26aOOoSzZtggED+KFhQwoeemjzSnPdsiuSeTSJHqaKFXl32inaWOqKMWNg9mzyHnqI404/fYsfqUCiSGZRDyRMKum+Y5Ytg7/8BU46Cc48M+poRGQ7lEDCpE2ltmuLBYNDhsD69bzWrx/FN90UbWAisl1KIGFSD2S7yhcMzh81Ch54gCV9+tDz0ku1YFCkDtAcSJjKE4jKmVQpFosx4aGHaNKtG98UFnLCM89owaBIHaEEEib1QJLS6YMPoLSU3qtXc9awYUoeInWEhrDCpASyfd9+yw/DhvFygwYcdM01WjAoUocogYRJk+jb9dFFF9Ho669pMno0111/vRYMitQhSiBhUg+kesuXs/vDD/PlCSdwxAUXAFowKFKXaA4kTJpEr97119OotJRd7rpri2YtGBSpG9QDCVPDhomHeiDbWrQI7rgDzj0XfvrTqKMRkRpQAgmbKvJW7oYboEGDxMpzEamTlEDCpgSyrcWL4f/+D847D9q2jToaEakhJZCwaVMpYKuSJTfeCA0a8PJxx1FcXBxtYCJSY0ogYcvP1yQ6P5YseeWRR2D8eJZ27cqpF1ygkiUidZjuwgqbhrCAH2/P/W/37hxZWkqPWbMoefxx3W0lUocpgYQtPx9Wrow6iowQ+9nPOH7jRsaVldFr8GAlD5E6TkNYYVMPZLNPLr2UeqWlrB08WCVLRLKAEkjYNIkOwKypUyksKWFlp0786dZbVbJEJAsklUDM7Cdm1iV4nWdm2p81WeqBALDxjjsodKfN8OGASpaIZIPtzoGY2bnAQKAlsC+wB3A7cGK4oWUJ3YUFpaWc+PbbcNxxcOSRm5tVskSkbkumB3IhcBywBsDdPwR2CTOobLB53UN+PmzYAKWlxOPx3Fz3MHEiLFkCl10WdSQikkLJJJAN7v5D+RszawB4eCFlh/J1D4uWLQNg1tNPU1RUlJvrHv71L9hnH+jVK+pIRCSFkkkgL5jZUCDPzH4JTACm1PbCZtbHzBaaWZmZdajmuK5m9r6ZLTKzKyu0721mrwftj5pZo9rGlErlY/x3PPAAABcMGEBJSUnuDdm89hq8+ipccgnUrx91NCKSQskkkCuBFcDbwHnAU8A1Kbj2O8BpwKyqDjCz+sAYoBtwANDXzA4IfjwcGOnu+wFfA+ekIKaUisViHN2lCwB/OPPM3EseALfeCgUFcPbZUUciIim23QTi7mXAA8A17n6Gu9/l7rUewnL399z9/e0cdiSwyN0XB8NojwCnmJkBnYHHguPuA06tbUypFo/HeWLmTACefPDBnLlldfP8z5dfwoQJcNZZxOfMyc35H5Estt0EYma9gAXAM8H7Q81scshxldsd+LTC+6VB287AN+6+aav2bZjZQDOba2ZzV6xYEWqwFcXjcYqKirjsuusAuPmaa3Jm3UP5/M//rr4aNm7k9cMOy935H5EslswQ1rUkegLfALj7AmDvZE5uZtPN7J1KHqfUNOAd5e53unsHd+/QunXrdF2WOXPmUFJSwuGdOwNwaLt2ObPuIRaLUfLwwzS4914++slP6PHnP+fm/I9IlkumFtZGd1+dGDXaLKkhLHfvUqOofvQZsGeF93sEbauA5mbWIOiFlLdnjCFDhiRefPJJ4nn1amJnnJEzX6KxH36AsjKKPv6YQcOG5cznFsklyfRAFprZb4D6ZtbezG4FXgk5rnJzgPbBHVeNgDOBycEcTBw4IzhuADApTTHtmMLCxPM330QaRrqtvP56lptx4NChqnslkqWSSSAXAQcCG4CHgNXAJbW9sJn1NrOlwDHAVDN7NmjfzcyeAgh6F4OBZ4H3gBJ3Xxic4grgUjNbRGJO5J7axhSKnXYCM1i9OupI0ubVhx+m5Wuv8X2/flx7ww2qeyWSpaodwgpuo53q7jHg6lRe2N0nAhMraf8c6F7h/VMkbh3e+rjFJOZmMlu9eonbWHOoB+J3343Vq0e7G24Atqx7paEskexRbQJx99JgoV+hu+fOn9Cp1rx57iSQ0lKOff996NoV9tprc7PqXolkn2Qm0dcCb5vZNGBzWVl3/2NoUWWb5s1zZwhr2jT47DMYPTrqSEQkZMkkkCeCR0WqhbUjCgtzpwdy772w887Qs2fUkYhIyJJJIM3dfYs/J83s4pDiyU7Nm8PHH0cdRfhWrYJJk2DQIGiUUaXJRCQEydyFNaCStrNSHEd2y5UeyEMPwQ8/qO6VSI6osgdiZn2B3wD7bFW6ZCfgq7ADyyq5Mol+771wxBFwyCFRRyIiaVDdENYrwDKgFTCiQvu3wFthBpV1mjeHNWugrCxxW282euMNWLAA/v3vqCMRkTSpMoG4+8fBQr/17v5CGmPKPoWF4A7ffvvjyvRsM348NG4MfftGHYmIpEm1fw67eylQZmZZ+q2XJs2bJ56z7FbezWXbN22CRx6Bnj2Jv/mmyraL5AitA0mH8gTyzTdbLK6r68rLtk+//HIO+fJL3j74YIqKiigpKYk6NBFJg5quA5EdkaUFFctLlHzQrRvtmzSh6+jRlEyYoBXnIjliuwnE3e9LRyBZLUuHsABiHTvyQ1kZ4zZu5Jw//1nJQySHJLMjYXsze8zM3jWzxeWPdASXNbK0BwLw7j//SaONG2kwYIDKtovkmGTuKR0HjAU2ATHgfhJ7pEuysrQHEo/H+fymm1i/yy6cc++9KtsukmOSSSB57j4DMHf/2N3/CpwcblhZJkt7IAvjcTqXltLknHOgXr0tyraLSPZLZhJ9g5nVAz40s8Ekto5tFm5YWaZRI8jLy7oEMrh168TiyP79N7epbLtI7kimB3Ix0BT4I3AE8Fsqr48l1cnGku6PPJIoW3LAAVFHIiIRSOYurPLxiLWAquTVVLYVVPz0U3jlFfjHP6KOREQiUl0xxSlUs++Hu/cKJaJslW09kPLFgkVF0cYhIpGprgdyc/BswF3AH8IPJ4s1b57YLyNblJQkKu/uu2/UkYhIRKorpri5gKKZrVVBxVoqLIT//S/qKFLjo49g9mwYPjzqSEQkQsnWFtcWtrWVTUNYGr4SEaqfA2lZ4W19M2tBYjgLAHfXplI7Ipsm0UtK4KijoF27qCMRkQhVNwcyj0TPozxpzK/wMwf2CSuorNS8eWK71/XroUmTqKOpuUWLYP58GDFi+8eKSFarbg5k73QGkvUqlnTfddcoI6md8uGrPn2ijUNEIpel+6tmoDpczmTzxlEAEybAMccQX7RIG0eJ5DglkHSpwwUVyzeOevXBB2HBAhYdcghFRUV07Ngx6tBEJELJ1MKSVKg4hFXHlBdJfKZHD44Bfv3oo5Q8/rhqXonkOCWQdCkfwqqDPRBIJJH9WrRg3rp1nDx4sJKHiOz4EJaZvRc8BocRUNaqwz0QgFdKStjzs8/4OhbTxlEiAtQggbj7z4HjgY9SH04Wq8OT6PF4nCm//z0AXW67TRtHiQiQ3Ja2FwWLCDdz91XuPjW8sLJQfj7Ur18nh7DmzJnDkPbtE2Xb999fG0eJCJBcD6QNMMfMSsysq5nZdn9DtmWWGMaqgz2QIWefTYu33oLTT9/cFovFGDJkSIRRiUjUtptA3P0aoD1wD3AWiZ0J/2FmKsO6owoL62QPhP/8J7HzYIUEIiKS1ByIuzvwRfDYBLQAHjMzrSTbEXW0B8ITTyTKth98cNSRiEgGSWYO5GIzmwcUAy8DB7n7IBLb29boT1Iz62NmC82szMw6VHNcVzN738wWmdmVFdrHm9lHZrYgeBxakzjSri4WVFyzBmbMgN69E8NwIiKBZNaBtAROc/ePKza6e5mZ9ajhdd8BTgPuqOoAM6sPjAF+CSwlMQ8z2d3fDQ75s7s/VsPrR6N5c/jww6ij2DFPPw0bN8Kpp0YdiYhkmGQSyGjYprz7t+6+0d3fq8lFy39vO/PxRwKL3H1xcOwjwCnAu9X9Ukari0NY//kP7LILHH101JGISIZJZg5kPrAC+AD4MHi9xMzmm9kRIca2O/BphfdLg7ZyN5jZW2Y20swaV3USMxtoZnPNbO6KFSvCijU5dW0SfcMGmDoVevVK3IIsIlJBMglkGtDd3Vu5+85AN+BJ4ALgtqp+ycymm9k7lTxOSUHcVwH7Ax1JDLFdUdWB7n6nu3dw9w6tW7dOwaVroXlz+PZb2LQp2jiSFY8n4tXwlYhUIpkhrKPd/dzyN+7+nJnd7O7nVfeXv7t3qWVsnwF7Vni/R9CGuy8L2jaY2Tjg8lpeKz3KV6OvWQMtW1Z/bCb4z38SCyBPPDHqSEQkAyXTA1lmZleY2U+CxxBgeTDJXRZibHOA9ma2t5k1As4EJgOYWdvg2YBTSUzKZ766VNK9rAwmTYJu3er2DooiEppkEshvSPz1/x9gIolewW+A+kBRTS5qZr3NbClwDDDVzJ4N2nczs6cA3H0TMBh4FngPKHH3hcEpHjSzt4G3gVbA32sSR9q1CCrCfFUHtpOfPRu++ELDVyJSpWqHsIJexmh371fFIYtqclF3n0giGW3d/jnQvcL7p4CnKjmuc02uG7lddkk8L18ebRzVKC4upmPHjsSefRYaNICTTyYejyfqYal0iYhUUG0CcffSYNiqkbv/kK6gslabNonnDE4g5bsPLsnLI79TJ+JvvEFRUREl5Xuhi4gEkplEXwy8bGaTge/KG939X6FFla3qQAKJxWJMuflm8s86iycPPJCzg+ShDaREZGvJJJD/BY96wE7hhpPl8vMTjwxOIABHB+tlLnzmGQYNG6bkISKV2m4Ccfe/AZhZU3dfF35IWa5Nm4xPIN/cfz+f1q/PgKFDGTt2LLFYTElERLaRTDHFY8zsXeC/wftDzKzKBYSyHRmeQF6aOJGd3n6bZv36cd1112n3QRGpUjK38Y4CTgJWAbj7m8AJIcaU3XbdNaMTyJqHH6Y+sPfFFwNo90ERqVIycyC4+6dbFT4sDSecHNCmDbz4YtRRVKn7pk2w++5w2GGb2zSEJSKVSaYH8qmZHQu4mTU0s8tJLOyTmmjTBlatysx6WOvXw7PPJoonau8PEdmOZBLI+cCFJCrhfgYcGryXmmjTBtwh6srAlZkxA9atg1NSUe9SRLJdMndhrQSqWokuO6riWpC2baONZWuTJ8NOO0GnTlFHIiJ1wHYTiJm1Bs4F2lU83t1/H15YWSxTFxOWlcGUKXDSSdC4yiLLIiKbJTOJPgl4EZiOJs9rrzyBfPFFtHFsbd48WLZMw1cikrRkEkhTd69ywybZQZnaA5k8ObHrYPfu2z9WRITkJtGfNDN9q6RKs2aQl5eZCeT44+vGRlcikhGSSSAXk0gi681sjZl9a2Zrwg4sa5ll3mr0JUvgrbcSt++KiCQpmbuwVEAx1TItgUyZknju2TPaOESkTkmmFpaZWX8zGxa839PMjgw/tCyWaQlk8mTYf39o3z7qSESkDklmCOs2ElvP/iZ4vxYYE1pEuSADEkhxcXGiQOLq1fD889CrF/F4nOLi4kjjEpG6I5kEcpS7XwisB3D3r4FGoUaV7dq0gZUroTS6u6LLdx5cOGIEbNrE/D32oKioiI4dO0YWk4jULcncxrsx2BvdYfPCwrJQo8p2bdokFu6tXPnjbb1pVl5l991u3WjXtCnd/vY3SiZMUNFEEUlaMj2QW4CJwC5mdgPwEvCPUKPKdhmyFiR2/PH0MGPCunWcd8EFSh4iskO2m0Dc/UFgCPBPYBlwqrtPCDuwrJYhCeSNW24hb/16Gvfpw9ixY7VplIjskGT3A/kvwY6EkgK77pp4jjCBxONxPrzmGg5u1Ii+48ax6+zZFBUVUVJSop6IiCQlmSEsSbUM6IHMmT2b3zZvTv1f/Qry87XzoIjsMCWQKBQUJCreRphAhnTrRt4XX2xRPDEWizFkyJDIYhKRukUJJM2Ki4uJP/98ohcSVOSNZP3FpEmJsipafS4iNaQEkmbl6y/WBAUV4/F4NOsvJk+Go46K7DZiEan7lEDSrHyu4dXFi1n25pvRTFx/9hnMnau9P0SkVpRAIhCLxWh14IGwfDmDBg1K/11PkycnnpVARKQWlEAiEI/HmfX+++xixu233Zb+9ReTJiUKJ+6/f3qvKyJZRQkkzcrnPE4+5xzqu/PEXXdRVFSUviSyejXMnJnofZil55oikpWUQNJszpw5lJSU8NPjjwfg+Pbt07v+4qmnYONG6N07PdcTkayV1Ep0SZ3N6yxefjnx/MknxLp3T988yBNPJFbCH310eq4nIllLPZCo7Ldf4vnDD9N3ze+/h6efhlNPhXr6Vy8itaNvkajssgvstFN6E8j06fDddxq+EpGUiCSBmFkfM1toZmVm1qGa4+41sy/N7J2t2lua2TQz+zB4bhF+1ClmlrgTKuQEsnnnQUgMXxUW8nzQLiJSG1H1QN4BTgNmbee48UDXStqvBGa4e3tgRvC+7klDAilf+f789OkwZQpfdOxIn379tPOgiNRaJAnE3d9z9/eTOG4W8FUlPzoFuC94fR9wauqiS6P27eHjj+GHH0K7RPnK91Gnnw6rVnHV66+rZLuIpERdnQNp4+7LgtdfAHWzoFP79omtbT/6KNTLxGIxrvzZz/ge2Ec7D4pIioSWQMxsupm9U8kjpfUz3N0J9muvIo6BZjbXzOauWLEilZeuvfbtE88hD2PFZ8xgz3nzWPLTn3LLPfdo50ERSYnQ1oG4e5ewzg0sN7O27r7MzNoCX1YTx53AnQAdOnSoMtFEIg0JJB6Pc9Ppp/NUWRm7/+UvlOy2m3YeFJGUqKtDWJOBAcHrAcCkCGOpuZ13hubNQ00gc+bM4Y4TT0xsYNWzp3YeFJGUscQIUJovatYbuBVoDXwDLHD3k8xsN+Bud+8eHPcw0AloBSwHrnX3e8xsZ6AE2Av4GChy98om27fQoUMHnzt3bgifqBaOPBIKC2HatHDOX1YGe+6ZuM7EieFcQ0SympnNc/dtllxEUsrE3ScC23ybufvnQPcK7/tW8furgBNDCzCd2rf/saxJGF55BT7/HIqKwruGiOSkujqElT3at4dPPoH168M5f0kJNGkCPXqEc34RyVlKIFFr3x7cYfHi1J+7tBQeewy6d0+UTRERSSElkKiFeSfWyy/DsmUavhKRUCiBRC3MBFJSAnl5cPLJqT+3iOQ8JZCotWiRuJ031Qlk06Yfh6+aNUvtuUVEUALJDGEUVXz2WVi+HPr3T+15RUQCSiCZIIwEct990KpVogciIhICJZBM0L49LF0K69al5nxffw2TJsFvfgONGqXmnCIiW1ECyQTlE+n/+1+tTrN586hHHkmUiB8wgHg8rs2jRCQUSiCZ4MADE8/z59fqNOWbR62+9VY46CDi33xDUVGRNo8SkVAogWSCAw9M3IlVyzLrsViMJ2++mcL33uOZNm0o+vWvVXVXREKjBJIJ6tWDTp0SCaSWxS2P+u9/KTXjrOnTGTRokJKHiIRGCSRTdO6cqIlVm5Immzax/u67mdGgAQOHDWPs2LHaPEpEQqMEkinKewq1+MJ/5/rrabJyJW3/8heuu+46SkpKKCoqUhIRkVAogWSK/feHXXetVQJpMX4863bbjYOuugpAm0eJSKiUQDJE8U03sfyAA2DmzM3zIDt0C+7rr7P7J5/Q9IoroH79zc2xWIwhQ4aEEbKI5DglkAzRsWNHbpw9G774At5/n3g8vmO34I4aBQUFcPbZocYpIlJOCSRDxGIxisaOBWDKpZdSVFSU/C24S5fChAnwhz9o3w8RSRslkAxyTL9+fFNQwPqnn96xW3D//e/EsNdFF4UboIhIBUogGST+/PM8u2ED3Zs25fbbbkvu7qmvv4Y774TevaFdu9BjFBEppwSSIcrnPA6+5BLy163jyRtvrPIW3M01rwCuvRZWr2bOSSep5pWIpJUSSIaYM2cOJSUl/PySS6BxY4586aUqb8Etr3k1++67YcwYlvbsSfehQ1XzSkTSyryWpTPqkg4dOvjcuXOjDmP7/vxnGDECFiyAgw+u9JD4zJk0PukkDmnYkEPz8rjzscdUtkREQmFm89y9w9bt6oFkoqFDoXlzuOKKKg+Jffklx27axJ++/56+F16o5CEiaacEkolatEgkkWeeSSws3NrixWwYPJg369en7dVXq+aViERCCSRTDR7M6ubNWTNoEJSVbW6eN2YMaw48kO+++oofbruNv/3976p5JSKRUALJVE2a8Nn551PwwQes3XdfuPxyPhw8mAMGD2ZDw4Z8MG4cHQcOBFTzSkSioUn0TFZWxvuXXsoXt93GcWVlNCgtZfUBB1D4/PPQunXU0YlIjtAkel1Urx4/GzWKGVdeSWFpKXf+/vcUzpun5CEiGUEJJMPF43HGjh3LZcOGcfXkycRffTXqkEREACWQjFa+Or2kpEQbRIlIxlECyWDlq9PL13hoslxEMokm0UVEpFqaRBcRkZRSAhERkRpRAhERkRqJJIGYWR8zW2hmZWa2zbhahePuNbMvzeydrdr/amafmdmC4NE9/KhFRKSiqHog7wCnAbO2c9x4oGsVPxvp7ocGj6dSGZyIiGxfgygu6u7vAZjZ9o6bZWbt0hGTiIjsmLo8BzLYzN4KhrlaVHWQmQ00s7lmNnfFihXpjE9EJKuFlkDMbLqZvVPJ45QUnH4ssC9wKLAMGFHVge5+p7t3cPcOrVVDSkQkZUIbwnL3LiGee3n5azO7C3gyrGuJiEjl6uQQlpm1rfC2N4lJeRERSaOobuPtbWZLgWOAqWb2bNC+m5k9VeG4h4FXgZ+Z2VIzOyf4UbGZvW1mbwEx4E9p/ggiIjlPtbBERKRaqoUlIiIplVM9EDNbAXwcvG0FrIwwnCjoM+eGXPvMufZ5If2f+Sfuvs1trDmVQCoys7mVdcmymT5zbsi1z5xrnxcy5zNrCEtERGpECURERGoklxPInVEHEAF95tyQa5851z4vZMhnztk5EBERqZ1c7oGIiEgtKIGIiEiN5FwCMbOuZva+mS0ysyujjicdqtrZMVuZ2Z5mFjezd4OdLy+OOqawmVkTM5ttZm8Gn/lvUceULmZW38zeMLOcKKpqZkuCUk4LzCzS0ho5NQdiZvWBD4BfAkuBOUBfd3830sBCZmYnAGuB+939F1HHE7ag2GZbd59vZjsB84BTs/nfsyV2Z8t397Vm1hB4CbjY3V+LOLTQmdmlQAegwN17RB1P2MxsCdDB3SNfPJlrPZAjgUXuvtjdfwAeAVKxP0lGc/dZwFdRx5Eu7r7M3ecHr78F3gN2jzaqcHnC2uBtw+CR9X8dmtkewMnA3VHHkotyLYHsDnxa4f1SsvyLJdcFWyIfBrwecSihC4ZyFgBfAtPcPes/MzAKGAKURRxHOjnwnJnNM7OBUQaSawlEcoiZNQMeBy5x9zVRxxM2dy9190OBPYAjzSyrhyvNrAfwpbvPizqWNDve3Q8HugEXBkPUkci1BPIZsGeF93sEbZJlgnmAx4EH3f2JqONJJ3f/BogDXSMOJWzHAb2COYFHgM5m9kC0IYXP3T8Lnr8EJpIYmo9EriWQOUB7M9vbzBoBZwKTI45JUiyYUL4HeM/d/xV1POlgZq3NrHnwOo/EjSL/jTSokLn7Ve6+h7u3I/H/8kx37x9xWKEys/zgxhDMLB/4FRHuyJpTCcTdNwGDgWdJTKyWuPvCaKMKXzU7O2ar44DfkviLdEHw6B51UCFrC8SDXTrnkJgDyYnbWnNMG+AlM3sTmA1Mdfdnogomp27jFRGR1MmpHoiIiKSOEoiIiNSIEoiIiNSIEoiIiNSIEoiIiNSIEoiIiNSIEoiIiNSIEohIDjCzW81svpl1jDoWyR5KICJZLih5sQtwHpD1+2VI+iiBiCTJzP5qZpcHr1+p5rjmZnZB+iKrnrt/R6LUyfPALdFGI9lECUSkBtz92Gp+3BzImARiZjsDTYFvgU0RhyNZRAlEpBpmdrWZfWBmLwE/q9C+NnjON7OpwV7k75jZr4EbgX2DIo43Bcf9J9gAaGH5JkBm1s7M3jOzu4L254JKupjZ78zsreC8/1fhuv2Dvc8XmNkdwTbN23MNcDOwEDgwRf9oRGgQdQAimcrMjiBRJvxQEv+vzCexv3pFXYHP3f3k4HcKSex++Itgc6dyv3f3r4IEMcfMHg/a2wN93f1cMysBTjezN0h86R/r7ivNrGVw7p8DvwaOc/eNZnYb0A+4v5rP0A44FrgUOJ5EAqly+E1kRyiBiFTt/wET3X0dgJlVtnfM28AIMxsOPOnuL5pZi0qO+6OZ9Q5e70kicXwBfOTuC4L2eUA7oAUwwd1XArh7+X72JwJHkEhAAHkktq+tzt+B69zdzew91AORFFICEakFd//AzA4HugN/N7MZbNUjMLNOQBfgGHdfZ2bPA02CH2+ocGgpiaRQFQPuc/erkonNzA4FTgOON7MxwTXfTuZ3RZKhORCRqs0CTjWzvGAXuJ5bH2BmuwHr3P0B4CbgcBKT1TtVOKwQ+DpIHvsDR2/nujOBPsHkN+VDWMAM4Awz26W83cx+EryeYWa7b3We4UAvd28X7Np3COqBSAqpByJSBXefb2aPAm+SGCqaU8lhBwE3mVkZsBEY5O6rzOxlM3sHeJrEfMb5wRDS+8Br27nuQjO7AXjBzEqBN4Cz3P1dM7sGeM7M6gXXu9DMPgX2A8qHujCzzkBTd59e4bzLzayZmbWsMCwmUmPakVCkjjOzX5CYpL806lgktyiBiIhIjWgOREREakQJREREakQJREREakQJREREakQJREREakQJREREakQJREREauT/AynF962/kyWlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(xdata, ydata, 'kx')\n",
    "x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)\n",
    "plt.plot(x, energy_surface.eval(x), 'r-')\n",
    "plt.xlabel(r'distance, $\\AA$')\n",
    "plt.ylabel('energy, Hartree')\n",
    "dist = max(ydata) - min(ydata)\n",
    "plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation of the molecular Vibronic Energy levels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Born-Oppeheimer approximation removes internuclear vibrations from the molecular Hamiltonian and the energy computed from quantum mechanical ground-state energy calculations using this approximation contain only the electronic energy. Since even at absolute zero internuclear vibrations still occur, a correction is required to obtain the true zero-temperature energy of a molecule. This correction is called the zero-point vibrational energy (ZPE), which is computed by summing the contribution from internuclear vibrational modes. Therefore, the next step in computing thermodynamic observables is determining the vibrational energy levels. This can be done by constructing the Hessian matrix based on computed single point energies close to the equilibrium bond length. The eigenvalues of the Hessian matrix can then be used to determine the vibrational energy levels and the zero-point vibrational energy \t\n",
    "\\begin{equation}\n",
    "{\\rm ZPE} = \\frac{1}{2}\\, \\sum_i ^M \\nu_i \\, ,\n",
    "\\end{equation}\n",
    "with $\\nu_i$ being the vibrational frequencies, $M = 3N − 6$ or $M = 3N − 5$ for non-linear or linear molecules, respectively, and $N$ is number of the particles.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we fit a \"full\" energy surface using a 1D spline potential and use it to evaluate molecular vibrational energy levels.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAELCAYAAAD3HtBMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4q0lEQVR4nO3deXhT1dbA4d+iTGUKQxkVrUIFEaQKRRSwDXIvoygOdcD7gRdFERwu4oCKA8hFigjiRZxbZyiKCoLK0CiiImUUUFFAHAFBBQRkatf3x0lCixRKydRkvc+TJ9mnJzkrJWT17H322qKqGGOMMceqTLgDMMYYUzpZAjHGGFMilkCMMcaUiCUQY4wxJWIJxBhjTImUDXcAoZSQkKCJiYnBP9CGDbBjB5x5ZvCPZYwxQbZkyZKtqlr70O0xlUASExNZvHhx8A80cCBkZ0MojmWMMUEmIt8fbrt1YQVDxYrw11/hjsIYY4LKEkgwxMfDnj3hjsIYY4LKEkgwVKwIeXmwf3+4IzHGmKCxBBIM8fHOvZ2FGGOimCWQYKhY0bm3cRBjTBSzBBIMlSo597t2hTcOY4wJIksgwVC5snNvCcQYE8UsgQSDJRBjTAywBBIMVao495ZAjDFRzBJIMNgZiDEmBlgCCQZfAtm5M7xxGGNMEFkCCQY7AzHGxABLIMFgCcQYEwMsgQSDDaIbc0wyMjLweDyFtnk8Hrp16/a37TfccAM33HBDwPctra+RkZFB2KhqzNxatWqlIZGfrxoXp3rPPaE5nokao0eP1pycnELb+vfvr/379y+0LScnR7t27fq3fYvaHumv0b9/f01ISPDvn5OTowkJCTp27Ni/ba9WrZq6XK6A71taX+PQ33EwAIv1MN+pYf9SD+XteBJIamaqZi7LVFXVfQf2aWpmqr684mVVVd21b5emZqbq5JWTVVV121/bNLVfnL45pLuqqm7ZtUVTM1N1+tfTVVV1458bNTUzVd/79j1VVf1h2w+ampmqc9bNUVXVdb+v09TMVP3wuw9VVfXrLV9ramaqfvLDJ6qqunLzSk3NTNVFPy1SVdVlG5dpamaqLtu4TFVVF/20SFMzU3Xl5pWqqvrJD59oamaqfr3la1VV/fC7DzU1M1XX/b5OVVXnrJujqZmp+sO2H1RV9b1v39PUzFTd+OdGVVWd/vV0Tc1M1S27tqiq6ptfvqmpmam67a9tqqo6eeVkTc1M1V37dqmq6ssrXtbUzFTdd2CfqqpmLsvU1MxU/+/ymcXP6AUvXuBvT1w0Ubu80sXfHv/ZeL3wtQv97TGfjNFLplzib4/6eJReMfUKf3v4h8O195u9/e1hOcO079t9/e2759yt10+/3t++/YPb9aZ3b/K3b33vVr31vVv97ZvevUlv/+B2f/v66dfr3XPu9rf7vt1Xh+UM87d7v9lbh3843N++YuoVOurjUarqJIQO/+ugYz4Z4//5uRPO1dP7nu7/j9/llS46cdFE/xdpuX7ldPCrg1XV+ZKI6xen8efGa05Oju47sE9bjmupVc6romPHjtVa9Wppy3EtdfLKyZqTk6M1G9TURiMaadU2VTUnJ0e37NqiLce11PjkeHVVq6Zvzpyiqc930IzXHtJTa9bUhx8crJWuK6vjXrhP9YcfNDtrnFa5rqzeeU9vbVOjhk55dpSmTkzRZ58cqm0rV9bkxHg9+7+n6aJ5L+mip57SNo2raNN7T9C2javo4okTddH7z+vZDzfW9qdW0dcHDdJzm1bWVg831q/fe0WXjB+vbZpU0Jr/LqPTx96r6vHopIwbtVa/OJ015j5dOnastj2zkja5q452qVdFlz76qE7PHqGthp+qXetU0cxrrtE2Z8drq+Gn6vwxI7RXtWp6443tNaFfnH6aMVJ1zhx9eNSVmtAvTp/tfZVeUq2aPvTfdG01/FS9pFo1zezdW89qV1HbDG+sOmeOLhszRpM7VNQWd9TXS6tV02UZGTo+a4CeP+J0vbRaNc3q3VtbuCtoxxHNdVlGhl5arZpeNihZT7yyrC7LyFCdPVsHjPynNryinGb17q2XVqum1z98gfad+E///t1vaaaJl5TTZRkZuiwjQ5N6lNeOtyX5j3fr07306ofb+4936sXl9F8Pn686e7Yuy8jQxEvKafdbmvn37zmitTbtUt5/vG4jztLhz/+f6uzZqrNn6xVPpOqozH6qc+aozpmjlzzRXsdk9XfaW7dqSVkCOc4E4hrl0n7v9FNV1V17d6lrlEtvnHGjqjoJwjXKpbfMukVVVb/f9r26horeMbiFqjoJwDXKpffNu09VVVdsWqGuUS4d8dEIVVVd+ONCdY1y6egFo1VV1fOdR12jXDr+s/Gqqjrrm1nqGuXSp3KfUlXVaV9OU9col2Yty1JV5wvcNcrlT2BZy7LUNcql076cpqqqT+U+pa5RLp31zSxVdb6gXaNc6vnOo6qqoxeMVtcoly78caGqqo74aIS6Rrl0xaYVqqp637z71DXK5U9Ad8y+Q12jXPr9tu9VVfWWWbeoa5TLn2BunHGjuka5dNdeJ6H0e6efuka5/L/Lf735L63xSA1/+4qpV2hCRoK/3WtyL607pq6/3f3V7tpgbAN/u/PLnbXhYw397Y5ZHTVxfKK/3eGFDtp4QmN/u+2zbbXpE0397VZPt9LmE5v728lPJWvyU8n+dvOJzbXV0wc/K02faKptn23rbzee0Fg7vNDB304cn6gdszr6zx4aPtZQO7/cWVVVx44dq9yOtp3gPD8nJ0flDtHmI5r7/3pMyEjQtAnna7OaNXVhZqZWe6iSnpxeVqf06qX/qVJFq95fXtMHN9Yn4uP1o5SztNLdaJ/r6qn26KHfnHOmVrkb/fellXR1XJx+nVRXXUNFB19YUX8T0RV1y6rrbvSejqiCrqiDuu5GR3Rw2gsbOO3R7Zy252SnPb6N057VyGk/1cppT2vqtLNaOu3JZzjtyWc47ayWTntaU6f9VCunPauR0x7fxml7Tnbao9s57YUNnPaIDk57RR2nfZ/baX9dy2nf0clpf1/Nad/S2WlviXfaN3Z32rvKOu1+PZ224tz+dTFa486D7SsuQxPuONjulY7WHXKw3f0qtMHgg+3OvdGGtx1sd/w/NPHWg+0O16KNbz7YbtsPbTrwYLvV9WjzAQfbyTc4N1+7+QBnH1+76UDnNXztxjc7x/C1E291YvC1G97mxOhrNxjsvAcF1ffe05IqKoHYGEgxJddLpv1J7QEoF1eO5HrJtDupHQCVylUiuV4y5zU8DwBXBRfJ2yrQdkc1AGpVqkVyvWTanNAGgDqV65BcL5nWDVoD0KBqA5LrJXN2/bMBOMl1Esn1kkmulwzAqTVOJbleMi3qtgAgqVYSyfWSaVa7GQBNEpqQXC+ZJglNAGhWuxnJ9ZJJqpUEQIu6LUiul8ypNU71v5fkesmc5DoJgLPrn01yvWQaVG0AQOsGrUmul0ydynUAaHNCG5LrJVOrUi0A2p7YluR6ybgquAA4r+F5JNdLplI5pwZYu5PakVwvmXJx5QBof1J7/3sB6HByB/97BTj/5PP9vwuA1JNT/b8rgLTENNqe2LZQ2/e79rXbNWxXqO37t/K1O5zcoVD7/JPPL3S81JNTC8WTlphWKN6Cbdd2F4kk+tvtGrbjZD2ZdevWkZ6eTuMKjUlLTMPj8TDuv/+lVfxp7HprGW9feCGLevSg0/46XLcB1tauTWKnTpy98jcuf3k+q3//nXOuvZaUdbu5d+0B0t96i8d27uTsH/bRNWct/ffvp/GSZZy1JY5OPwr8/DMnlIun2baKJHy9G01Kol6b9iRrHdqd0oFvW7XinT8PcMruarQ993J44AHWJJ9LtU0gdZNh3Dga3PEQyVWT2NOwDdcBX7XuSLKrCcn97oFXX2VNSifKb4LNzdLgjTdIGjaOBnm1eX4LvJqeTpNh40lOaE6TBybw8pVX8vQWqJ9fh6T/PgUffECLoeOol1+H4dsg65prSL5rLMl1W/J5x4txA7+e0Ynkui1p8PwUlo0bx8ztlaindRkSV5Wl48bRZtB/aVy5EbeUrcrz//d/5PxRicaVG/HdQ6O40OViT6PzidsSx5pRY+Djjzkh9RritsTxUu9/0cPl4sT2V9O4SmN6uFw816cPy7fEk1TtNPj4Y5Y+/jirNleknpxAd5eLpRMmkHrxLTSp0ZzuLhfP9e3Lus0VaVqrJUsnTKC7y0XlhLM5sLkcSydMgAULOO2MbuzbVI7n+valu8tFkyadad+8m3//GlVbsGuTs//SCRP4fWN5EuKbHjxem8s5PTHVf7wdG8tzeqOOsGABSydMYNemctSo2sK//+n1z2Hjxgr+451etzVpna6DBQtgwQLOO60jaV1ugI8/ho8/pm2jVNK63+S02xz8PxUwh8sq0XoL2RiIqmrLlqoXXnjU3UzkOtx4xGH76+fN0xY1aujyUaP02wED9OUKFfS7k07SX8qU8f8lWOiWkKDaooVqp066vHlzHQc6Ny1NdeJE1ddf1+WPPKKdXS59/MYbtXmNGvrRjBmaM2eOJiQk6LBhww7bZ16c7ceyb6hf49D+fBsDsTGQiLuFNIG0a6fasWPojmcCrqgvtY+nTdMVI0bo+Ph4XXvKKbpFpFCC+LNSJZ0PuvTMM1WHD9fVQ4fqhS6Xjr35Zq1fq1axvlxL0xdYIF6jf//+QRnMj/SLBwLxGqNHj9ZgswRynAnkmAfR/1ND3+yZpKo2iF4aBtHbjWjn/w/rG0T3nW1UvayC9hzYVJ+vWFF3NmyofS9Ch7nRA2XK6GLQLv1q6PCM7qrz5+v8t97S8leX1473d/R/iZb/V3m94cUbVNX5D1++T3m9aNRF/i/XLq900VtfudX/RXr2+LP1mcXP+OOrP7S+pt7q/P58n72hrw/Vrl276qw5swp99qbPnq41B9fUB7MfVNWDn73ON3fW/v37F/rs5eTkaNpFadpyXMtCn72W41pqm/Q2mpOTU+iz179/f73spssKffaemf6M1hxcU5+Z7sTr++w9P+N57dq1q054a0Khz16PQT20/tD6hT57Lce11KGjhsbsZy9YF3Coql4y5ZJCF3CUVFEJxMZAgiUuzlYkLEXq1qtLenq6/9r7X9asYfaFFzJq8WKu2LuXZl9/Te+8PCq3aAGtW7Gxc2cSa9TgnWHD8OzezXe16+A5cIBLrr+edu3a8Y9//IPs7GyGDRtGkyZNaNy4MQBut5tWrVrxzZpvyM7Oxu12A3DaaaeRnZ1No0aNqFGjRqHYTjvtNPr27VtoW7NmzZg1axapqamFtp/f4XxatGhBi+YtCm0fOHAgTz/9dKFtbrebl158ierVqxfaXr16dTJGZ/hj83n66ad54IEHCm1LSUmhRYsWpKSkFNrepk0bZs2aRatWrQptH3L7EE477bS/HW/AgAGYUuhwWSVabyHtwrrmGtXExNAdzxy3j2bM0JurVNF1iYma5+2O2nnSSfp4fLw+16ePv/upqK6torphQtHFYEwwYWcgIVa5ss1Ej0CHm/H8eVYWS887j/OvuooJO3eSv2EDH55/Pp9nZZG4ezctZs6kX1YWr06dSnp6OpMnTy509uB2u/1nD4f+1e52u7nzzjtD9v6MCaWy4Q4galWpYgkkAqWkpJCenu4kgHLl+P0//+GcxYvJK1eOjR078u+FC0m5+WYmPfUUl3z66WETRW5u7mETxaHbjIl2lkCCpXJl2L0b8vOhjJ3oRQq3280H993Hrs6dYf9+Doiwvl8/funenV79+5P91ltOMujYkfT0dK688sq/Pd8ShTEO+2YrprSsNLKWZwGwP28/aVlpvPLFKwDs3r+btKw0pqyaAsD2PdtJK/sK004H/vqLrbu3kpaVxow1MwDYtHMTaVlpvL/2fQB+3P4jaVlpzF0/F4D1f6wnLSuNjzZ8BMCarWtIy0rj0x8/BWDVr6tIy0oj9+dcAJZvWk5aVhrLNy0HIPfnXNKy0lj16yoAPv3xU9Ky0lizdQ0AH234iLSsNNb/sR6AuevnkpaVxo/bfwTg/bXvk5aVxqadmwCYsWYGaVlpbN29FYBpX00jLSuN7Xu2AzBl1RTSstLYvX83AK988QppWWnsz9sPQNbyLNKy0vy/y2eXPEunlzr520/mPknXV7v6248vfJyer/f0tx/99FEuzb7U335kwSNc+cbBL/YRH43gmmnX+Nv3e+6n9cOt/V1VQ+cOpf+M/nz22musadGC12bdxkvd4Dbg6bvuYsJlVbh79Vj/2cbAmQOZuX+m/2yj/4z+DJ071P/6175zLfd77ve3r5l2DSM+GuFvX/nGlTyy4BF/+9LsS3n000f97Z6v9+TxhY/7211f7cqTuU/6251e6sSzS571t4/5s5eVxrSvpgHYZy8Mn71r37nW3/Z99nyGzB7CwJkD/e3b3r+N296/zd8eOHMgQ2YP8bcD/dkLNEsgweI767BFpcIiISHh4FVVeQf40+PhrN69abx2LdvOOIPX48pQbdgwJjz3HD/+9COtW7e28QtjjtXhRtaDfQNqAnOAb733NYrYbzSwynu7osD2LOA7YLn3llyc44b0KqzMTFVQXbcudMc0heTk5Og/XC7dnJCgCrqxUyddMGXKYa+gCsVsXmNKKyLsKqy7gXmqmgTM87YLEZHuwNlAMnAOMEREqhXY5Q5VTfbelgc/5GNki0qF1/79uOfM4f0dO9i7dSsvXXUV9ebM4ZMNG4ocGDfGHJtwDaJfBKR5H78IfAjcdcg+zYD5qnoAOCAiXwBdgOwQxXh8bFGpkMnIyCAlJeVgF9Qvv7CtSxeqr1zJ6xUqsOGWWxifmUlDj+ewXVI2MG5MyYTrDKSuqm70Pt4E1D3MPiuALiJSSUQSADfQsMDPR4rIFyIyTkQqFHUgEekvIotFZPGWLVsC9gaOys5AQsZ3aa7H44GPPmJv8+aUW7mS6+LjafDee9ybkUF2dnahmebGmOMXtAQiInNFZNVhbhcV3M/bv6aHPl9VZwOzgE+B14HPgDzvj4cCTYEUnPGUQ89eCr7OM6raWlVb165dOyDvrVh8CcQG0YPO1w31Vs+e5Lnd/LhjB2Muu4zeM2daV5UxQRS0LixV7VTUz0Rks4jUV9WNIlIf+LWI1xgJjPQ+5zXgG+9239nLXhHJBIYc7vlhZWcgoaOK+5NPcO/cyWwg9/bbeXD06L/tZl1VxgRWuLqwpgN9vI/7AO8cuoOIxIlILe/jM4Ezgdnedn3vvQAX41ylFVksgYRGXh4MGADDhjGlQgU+u+cexr/wgnVVGRMC4UogjwD/EJFvgU7eNiLSWkSe8+5TDvhYRL4EngGu8Q6oA7wqIiuBlUAC8HCwAz7myVwfXOVMJNy1yyZzBWsyV34+XH89I756mtaXx1Fn1iweGDmSrhld6fp0179NJPQJ92Qum0gYBZ89r1ifSBiWq7BU9TfggsNsXwxc5328B+dKrMM9v2NQAwyEuDjn3s5AAsZ3tRXlcGrl3ngjZGayoldNarjPxt3R+VgkJiaSlpZ22JpVxpjAEWcMOza0bt1aFy9eHLoDlisHQ4bAqFGhO2YU83g8TiHEKVNwv/kmPPkk4+LjSX73XX/yMMYEnogsUdXWh263YorBZCXdA8p3JdXCHj1w797NE5Y8jAkrSyDBZAkk4Ny//IJ7925eBLbcfrslD2PCyIopBpMlkMD6+GPyr72WBeXKseGee5j01FN2tZUxYWQJJJhsUanA+e479vfowbr8fHTqVB4YOdJmlxsTZpZAgsnOQAJjzx64/HLy9u9na1YWHS5yihnY7HJjwssSSDBVrmylTEqo0Nrlt98OS5bwzT338PEvvxTaz9bsMCZ8LIEUU4kmc9XfZhMJKdlkrrfj3yY9PZ3Vw4bx6PIn6Tr4BC54/HFSUlJK/WQum0gY2Z89m0hYfJZAgqliRevCKqGEhARmjBvHySNH8ofLhefXLYXW8TDGhJ9NJAymm26CqVMhlGXko0VeHqSm8tfixZy2dy/XDhvG8OHDwx2VMTGpqImEdgYSTDaIXnKPPw6ffMKQcuW4dtgwJk2aZFdbGRNhLIEEU+XK8Ndfzl/Tpvi+/pq8u+/mvfLlueyddxg+fLhdsmtMBLIEEky+ku67d4c3jtLkwAHo04d95cpR7bXX/DPN7ZJdYyKPlTIJpoJrglStGt5YSouJE2HRIuJfe412l15a6Ee2IJQxkcXOQILJFpU6Nhs3wv33Q+fOcOWVR9/fGBNWlkCCqUoV594SSJEKTRi8807Ys4eFvXuTMWZMeAMzxhyVJZBgsjOQo0pJSSE9PZ2l48fDK6+w4fLLuXDwYGfhKGNMRLMxkGDyJRArZ1Ikt9vN1Ndeo2LXrmxzuTj//ffJnjrVxjqMKQUsgQSTnYEUS9o330BeHr22b6fvsGGWPIwpJawLq5hKVI/o908A2Lpjk9UjKqoe0Z9/MmHqHbS9Rmhx331MmjSJQa8Miup6RFYLK0I+e1gtrONlCSSYKlZ07m0eSJG+u/lmyu7+i7jTT2f4iBFkZ2eTlZnF1q1bwx2aMeYorBZWMG3bBjVqwGOPwX/+E7rjlhabN7PvpJPY1rYtdT76yL/Z4/GQm5trZdqNiRBF1cKyMZBgskH0IxsxgvJ5edR59tlCm23CoDGlg3VhBVO5cs7NBtH/bu1aePppuP56OO20cEdjjCkBSyDBZhV5D2/kSChb1pl5bowplSyBBJslkL9bvx5efhluuAHq1w93NMaYErIEEmxVqlgC4ZCSJY88AmXL8km7dmRkZIQ3MGNMiVkCCbbKlW0QnYMlSz6dPBmysvipSxcuvukmK1liTClmV2EFm3VhAQfX8/i6Wzfa5OXRY/58st980662MqYUswQSbJUrg02KA8DdpAnt9+8nMz+fnoMGWfIwppSzLqxgszMQvx8GD6ZMXh47Bw2yNc6NiQKWQILNBtEBmD9zJq7sbLampfGfJ56wNc6NiQLFSiAicrKIdPI+jheRmFuftcQF7SpXZmvenzFf0O5mTx9cqtQdPZpnlzzLyB9H+tc4j7WCdlZM0Yop+kR9MUURuR54A3jau+lE4O2gRRRtrAsL8vOpsWsntGsHbdr4N7vdbqt3ZUwpdtRiiiKyHGgDfK6qZ3m3rVTVFsEPL7BCWUwxIyODlJQU3PPnw4MPwoEDeObPj80igW+8AZdfDtOmQa9e4Y7GGHOMiiqmWJwurL2quq/AC5UFYqeEbwn55j2s3bgRgPnvvUd6enpsznt47DE49VTo2fPo+xpjSo3iJJCPROQeIF5E/gFMBWYc74FF5HIRWS0i+SLyt8xWYL8uIrJGRNaKyN0Ftp8iIp97t08RkfLHG1Mg+eY9PP2K01d9U58+ZGdnx96lqwsXwmefwW23QVxcuKMxxgRQcRLI3cAWYCVwAzALuC8Ax14FXALML2oHEYkDJgJdgWbAVSLSzPvj0cA4VW0M/AH0C0BMAeV2u2nbyRmwu+7KK2MveQA88QRUqwbXXnv0fY0xpcpRE4iq5gOvAPep6mWq+qwGYBUqVf1KVdccZbc2wFpVXe/tRpsMXCQiAnTEGdwHeBG4+HhjCjSPx8O0nBwA3n311Zi5ZNVf9+rXX2HqVOjbF09urtW9MibKFOcqrJ7AcuB9bztZRKYHOS6fE4AfC7R/8m6rBWxT1QOHbP8bEekvIotFZPGWLVuCGmxBHo+H9PR0bh8+HIBH77svZuY9+MZ/1t17L+zfz+dnnRW74z/GRLHidGE9gHMmsA1AVZcDpxTnxUVkroisOsztopIGfKxU9RlVba2qrWvXrh2qw5Kbm0t2djZnd+wIQHJion/eQ7Rzu91kv/46ZV94ge9OPpked9wRm+M/xkS54tTC2q+q251eI79idWGpaqej73VEPwMNC7RP9G77DaguImW9ZyG+7RHDf6nuDz8499u3477sspj5EnXv2wf5+aR//z0Dhg2LmfdtTCwpzhnIahG5GogTkSQReQL4NMhx+eQCSd4rrsoDVwLTvWMwHuAy7359gHdCFNOxcbmc+23bwhpGqG0dMYLNIpxxzz1W98qYKFWcBHIzcAawF3gN2A7cdrwHFpFeIvITcC4wU0Q+8G5vICKzALxnF4OAD4CvgGxVXe19ibuAwSKyFmdM5PnjjSkoqlYFEdi+PdyRhMxnr79OzYUL+at3bx4YOdLqXhkTpY7YheW9jHamqrqBewN5YFV9C3jrMNt/AboVaM/CuXT40P3W44zNRLYyZZzLWGPoDESfew4pU4bEkSOBg3NicnNzrSvLmChyxASiqnneiX4uVY2dP6EDrXr12EkgeXmct2YNdOkCJ53k3+x2uy15GBNlitOFtRNYKSLPi8gE3y3YgUWa46qIWqcKaXVmxUZF1DlzeKXWz6SlbbCKqF5Wjdeq8fqE+7MXaMW5Cmua91aQ1cI6FtWqwYGfwh1FaLzwgrMGSkJCuCMxxgRZcarx3qqqjx9tW2kQymq8hVx0EXz/PSxfHvpjh9Jvv0GDBjBgAIwfH+5ojDEBcjzVePscZlvf444olrhcsTEG8tprsG+f1b0yJkYU2YUlIlcBVwOnHlK6pCrwe7ADiyqxMoj+wgvQqhW0bBnuSIwxIXCkMZBPgY1AAjC2wPY/gS+CGVTUqV4dduyA/Hznst5otGyZ00X3v/+FOxJjTIgUmUBU9XvvRL89qvpRCGOKPi4XqMKffx6cmR5tsrKgQgW46qpwR2KMCZEj/jmsqnlAvohE6bdeiFSv7txH2Wx0f9n2Awdg8mS48EI8K1ZY2XZjYkRxLuP1zQOZA+zybVTVW4IWVbTxJZBt2wpNrivtfGXb5w4ZQstff2XlmWeSnp5OdnZ2uEMzxoRASeeBmGMRpQUVfSVKvunalaSKFeny+ONkT51qM86NiRFHTSCq+mIoAolqUdqFBeBOSWFffj6Z+/fT7447LHkYE0OKsyJhkoi8ISJfish63y0UwUWNKD0DAfhy1CjK799P2T59rGy7MTGmONeUZgKTgAOAG3gJZ410U1xRegbi8Xj4ZcwY9tSpQ78XXrCy7cbEmOIkkHhVnYdT9uR7VX0Q6B7csKJMlJ6BrPZ46JiXR8V+/aBMmUJl240x0a84g+h7RaQM8K2IDMJZOrZKcMOKMuXLQ3x81CWQQbVrO5MjrzlYjdTKthsTO4pzBnIrUAm4BWgF/IvD18eKasddUvtfeczY60zgj5qS2m+8Ai1bMiV/pZXUtnLugJVzj/TPXqAV5yosX3/ETsCq5JVU2TjYtevo+5UWe/fC54vggf+GOxJjTJgUWc5dRGZwhHU/VLVnUT+LVGEr5w5w7rnO+uizZ4fn+IE2diwMGQJr10KjRuGOxhgTREWVcz/SGYjvvEeAZ4HrghFYzKhe3VkvI1pkZzuVdy15GBOzjlRM0V9AUUR2WkHF4+Rywbp14Y4iML77DhYtgtGjwx2JMSaMiltb3JawPV7Vq0fPPBBfrav09PDGYYwJqyMtKFWzQDNORGrgdGcBoKq2qNSxiKZVCbOz4ZxzIDEx3JEYY8LoSGMgS3DOPHxJY2mBnylwarCCikrVqzvLve7ZAxUrhjuaklu7FpYudQbRjTEx7UhjIKeEMpCoV7Cke7164Yzk+Pi6ry6/PLxxGGPCLkrXV41ApbiciX/hKICpU+Hcc/GsXWsLRxkT4yyBhEopLqjoWzjqs1dfheXLWduyJenp6aSkpIQ7NGNMGBWnFpYJhIJdWKWMr0ji+z16cC5wxZQpZL/5ptW8MibGWQIJFV8XVik8AwEniTSuUYMlu3fTfdAgSx7GmGPvwhKRr7y3QcEIKGqV4jMQgE+zs2n488/84XbbwlHGGKAECURVTwfaA98FPpwoVooH0T0eDzP+/W8AOj35pC0cZYwBirek7c3eSYR+qvqbqs4MXlhRqHJliIsrlV1Yubm53JmUBM2aQdOmtnCUMQYo3hlIXSBXRLJFpIuIyFGfEYWOe02GF93MSK4E27aVujUZTnc3oFfycrZe6qybMO2raTz0/UPccMsNgK3JYOuB2HogPpH+2Qu0oyYQVb0PSAKeB/rirEz4XxGxMqzHqnKlUnkGwsKFzn2PHuGNwxgTUYpcD+RvO4q0xFlQqgvgAdoCc1T1zuCFF1hhXQ8EnPLn9evDu++GL4aS6NoVvv3WucXmCagxMa2o9UCKMwZyq4gsATKAT4AWqjoAZ3nbS4/45KJf83IRWS0i+SLyt6AK7NdFRNaIyFoRubvA9iwR+U5ElntvySWJI+RKY0HFHTtg3jzo1cuShzGmkOLMA6kJXKKq3xfcqKr5IlLSPo1VwCXA00XtICJxwETgH8BPOOMw01X1S+8ud6jqGyU8fnhUr+78FV+avPce7N8PF18c7kiMMRGmOAnkcfhbefc/VXW/qn5VkoP6nneU8fg2wFpVXe/ddzJwEfDlkZ4U0apXL31nIG+/DXXqQNu24Y7EGBNhinMV1lJgC/AN8K338QYRWSoirYIY2wnAjwXaP3m3+YwUkS9EZJyIVCjqRUSkv4gsFpHFW7ZsCVasxeNyla5B9L17YeZM6NnTuQTZGGMKKE4CmQN0U9UEVa0FdAXeBW4CnizqSSIyV0RWHeZ2UQDiHgo0BVJwutjuKmpHVX1GVVurauvatWsH4NDHoXp1+PNPOHAgvHEUl8fjxGvdV8aYwyhOF1ZbVb3e11DV2SLyqKrecKS//FW1U1E/K6afgYYF2id6t6GqG73b9opIJjCE0sA3G33HDqhZ88j7RoK333YmQF5wQbgjMcZEoOKcgWwUkbtE5GTv7U5gs3eQOz+IseUCSSJyioiUB64EpgOISH3vvQAX4wzKR77SVNI9Px/eece5hLc0r6BojAma4iSQq3H++n8beAvnrOBqIA5IL8lBRaSXiPwEnAvMFJEPvNsbiMgsAFU9AAwCPgC+ArJVdbX3JV4VkZXASiABeLgkcYRcDW9FmN9LwXLyixbBpk3WfWWMKdIRu7C8ZxmPq2rvInZZW5KDqupbOMno0O2/AN0KtGcBsw6zX8eSHDfs6tRx7jdvDm8cR5CRkUFKSgruDz6AsmWhe3c8Ho9TD+vOUjNn1BgTAkdMIKqa5+22Kq+q+0IVVNSqW9e5j+AE4lt9cEN8PJXT0vAsW0Z6ejrZvrXQjTHGqziD6OuBT0RkOrDLt1FVHwtaVNGqFCQQt9vNjEcfpXLfvrx7xhlc600etoCUMeZQxUkg67y3MkDV4IYT5SpXdm4RnEAA2nrnywx8/30GDBtmycMYc1hHTSCq+hCAiFRS1d3BDynK1a0b8Qlk20sv8WNcHH3uuYdJkybhdrstiRhj/qY4xRTPFZEvga+97ZYiUuQEQnMUEZ5AFrz1FlVXrqRK794MHz7cVh80xhSpOJfxjgc6A78BqOoK4PwgxhSRAraoT716bNr2U8Qu6rPj9deZ0xiuTVnNpp2bcLvd3DbpNq5bcJ0t6uNlC0rZglI+pe2zF2jFWhNdVX88ZFNeEGKJDXXrwtat4Y6iSN0OHIBataBKFf+2M1ucScOGDY/wLGNMLDrqglIi8gbwGPA/4BzgVqC1ql55xCdGoLAvKAXw4IMwfDjs2+fMs4gke/Y4yaNPH3jSeimNMY4SLygF3AgMxKmE+zOQ7G2bkqhbF1Qh3JWBD2fePNi9Gy4KRL1LY0y0K85VWFuBomaim2NVcC5I/frhjeVQ06dD1aqQlhbuSIwxpcBRE4iI1AauBxIL7q+q/w5eWFEsUicT5ufDjBnQuTNUKLLIsjHG+BWnE/4d4GNgLjZ4fvx8CWTTpvDGcaglS2DjRuu+MsYUW3ESSCVVLXLBJnOMIvUMZPp0Z9XBbt2Ovq8xxlC8QfR3RcS+VQKlShWIj4/MBNK+felY6MoYExGKk0BuxUkie0Rkh4j8KSI7gh1Y1BKJvNnoGzbAF184a58bY0wxFecqLCugGGiRlkBmOLOUufDC8MZhjClVilMLS0TkGhEZ5m03FJE2wQ8tikVaApk+HZo2haSkcEdijClFitOF9STO0rNXe9s7gYlBiygWREACycjIcAokbt8OH34IPXvi8XjIyMgIa1zGmNKjOAnkHFUdCOwBUNU/gPJBjSra+eph5YXvqmjfyoOrx46FAwdYeuKJpKenk5KSEraYjDGlS3Eu493vXRtdwT+xMD+oUUW7unWdiXtbtx68rDfE3G432dnZfNm1K4mVKtH1oYfInjrV1v0wxhRbcc5AJgBvAXVEZCSwAPhvUKOKdhEyF8Tdvj09RJi6ezc33HSTJQ9jzDE5agJR1VeBO4FRwEbgYlWdGuzAolqEJJBlEyYQv2cPFS6/nEmTJtmiUcaYY1KseuKq+jXeFQlNANSr59yHMYF4PB6+ve8+zixfnqsyM6m3aBHp6elkZ2fbmYgxpliKtaCUCfCqcFUgrS+8/7Oz6ls4VoW77uN+tG9Yhbh//pO5mz/joe8f4n8v/Y/c3FxbFc5WJLQVCb2i7bMXaJZAwqFqVWdG+rZtYQuhYa1aVNy6tVDxxPPOO48777wzbDEZY0qXo65IGE0iYUXCjIwMUlJScPftC6mp8NJLeDwecnNzQ/vlPWIEPPCAU4E3TFeCGWNKh+NZkdAEkG/+xQ5vQUWPxxOe+RfTp8M551jyMMaUmCWQEPPNv/hs/Xo2rlgRnoHrn3+GxYtt7Q9jzHGxBBIGbrebhDPOgM2bGTBgQOivepo+3bm3BGKMOQ6WQMLA4/Ewf80a6ojw1JNPhn7+xTvvOIUTmzYN7XGNMVHFEkiI+cY8uvfrR5wq0559lvT09NAlke3bISfHOfsQCc0xjTFRyRJIiOXm5pKdnc1p7dsD0D4piezsbHJzc0MTwKxZsH8/9OoVmuMZY6JWsWaim8DxX6r7ySfO/Q8/4O7WLXTjINOmOTPh27YNzfGMMVHLzkDCpXFj5/7bb0N3zL/+gvfeg4svhjL2T2+MOT72LRIudeo4M9JDmUDmzoVdu6z7yhgTEGFJICJyuYisFpF8Efnb7MYC+70gIr+KyKpDttcUkTki8q33vkbwow4wEedKqCAnEP/Kg+B0X7lcfOjdbowxxyNcZyCrgEuA+UfZLwvocpjtdwPzVDUJmOdtlz4hSCC+me8fzp0LM2awKSWFy3v3tpUHjTHHLSwJRFW/UtU1xdhvPvD7YX50EfCi9/GLwMWBiy6EkpLg++9h376gHcI38338pZfCb78x9PPPrWS7MSYgSusYSF1V3eh9vAkonQWdkpKcpW2/+y6oh3G73dzdpAl/AafayoPGmAAJWgIRkbkisuowt4DWz1CnnHCRJYVFpL+ILBaRxVu2bAnkoY9fUpJzH+RuLM+8eTRcsoQNp53GhOeft5UHjTEBEbR5IKra6eh7ldhmEamvqhtFpD7w6xHieAZ4Bpxy7kGM6diFIIF4PB7GXHops/LzOeH++8lu0MBWHjTGBERp7cKaDvTxPu4DvBPsAwZlVbg/cqF6dX5ctyxoq8Ll5uZy/YXJpF0rrGnXFLfbzb3P3Mt1C66zVeG8bEVCW5HQJ9o/e4EWrst4e4nIT8C5wEwR+cC7vYGIzCqw3+vAZ0ATEflJRPp5f/QI8A8R+Rbo5G2XPr5LeYM4BnLnkCHUXb0aatWEKlUAOOuss2jYsGHQjmmMiQ22ImG49e7tlDXZsCE4r79gAXToAK+9BlddFZxjGGOimq1IGKmSkuCHH2DPnuC8fnY2VKwIPXoE5/WNMTHLEki4JSWBKqxfH/jXzsuDN96Abt2csinGGBNAlkDCLZhXYn3yCWzcCOnpgX9tY0zMswQSbsFMINnZEB8P3bsH/rWNMTHPEki41agBtWoFPoEcOHCw+8p79ZUxxgSSJZBIEIyiih98AJs3wzXXHH1fY4wpAUsgkSAYCeTFFyEhwTkDMcaYILAEEgmSkuCnn2D37sC83h9/wDvvwNVXQ/nygXlNY4w5hCWQSOAbSF+37rhexr941OTJTon4Pn3weDy2eJQxJigsgUSCM85w7pcuPa6X8S0etf2JJ6BFCzzbtpGenm6LRxljgsISSCQ44wznSqzjLLPudrt599FHcX31Fe/XrUv6FVdY1V1jTNBYAokEZcpAWpqTQI6zNtk5X39Nngh9585lwIABljyMMUFjCSRSdOzo1MQ6npImBw6w57nnmFe2LP2HDWPSpEm2eJQxJmgsgUQK35nCcXzhrxoxgopbt1L//vsZPnw42dnZpKenWxIxxgSFJZBI0bQp1Kt3XAmkRlYWuxs0oMVQZwEat9tNdnY2ubm5gYrSGGP8LIFEiIwxY9jcrBnk5PjHQY7pEtzPP+eEH36g0l13QVycf7Pb7ebOO+8MRsjGmBhnCSRCpKSk8MiiRbBpE6xZg8fjObZLcMePh2rV4Nprj7qrMcYEgiWQYgr2utQPff8Qp4x11kp+duiN/PP1f3LvM/fidruPvi71ig9IqziZ5f17QtWqti51hK9LbWui22fPx9ZENwHT7IIL2FatGns/+ogTGpzAWWedVbwnTp7s3NuStcaYELI10SOIx+Ph165d6REXxynx8UyZOvXo8zj++AMaNXIuA37jjdAEaoyJKbYmeoTzjXmcedttVN69m3cfeaTIS3D9Na8AHngAtm8nt3Nnq3lljAkpSyARIjc3l+zsbE6/7TaoUIE2CxYUeQmur+bVoueeg4kT+enCC+l2zz1W88oYE1LWhRWJ7rgDxo6F5cvhzDMPu4snJ4cKnTvTslw5kuPjeeaNN6xsiTEmKKwLqzS55x6oXh3uuqvIXdy//sp5Bw7wn7/+4qqBAy15GGNCzhJIJKpRw0ki77/vTCw81Pr17B00iBVxcdS/916reWWMCQtLIJFq0CC2V6/OjgEDID/fv3nJxInsOOMMdv3+O/uefJKHHn7Yal4ZY8LCEkikqliRn2+8kWrffMPORo1gyBC+HTSIZoMGsbdcOb7JzCSlvzNByWpeGWPCwQbRI1l+PmsGD2bTk0/SLj+fsnl5bG/WDNeHH0Lt2uGOzhgTI2wQvTQqU4Ym48cz7+67ceXl8cy//41ryRJLHsaYiGAJJMJ5PB4mTZrE7cOGce/06Xg++yzcIRljDGAJJKL5ZqdnZ2fbAlHGmIhjCSSC+Wan++Z42GC5MSaS2CC6McaYI7JBdGOMMQFlCcQYY0yJWAIxxhhTImFJICJyuYisFpF8Eflbv1qB/V4QkV9FZNUh2x8UkZ9FZLn31i34URtjjCkoXGcgq4BLgPlH2S8L6FLEz8aparL3NiuQwRljjDm6suE4qKp+BSAiR9tvvogkhiImY4wxx6Y0j4EMEpEvvN1cNYraSUT6i8hiEVm8ZcuWUMZnjDFRLWgJRETmisiqw9wuCsDLTwIaAcnARmBsUTuq6jOq2lpVW9e2GlLGGBMwQevCUtVOQXztzb7HIvIs8G6wjmWMMebwSmUXlojUL9DshTMob4wxJoTCdRlvLxH5CTgXmCkiH3i3NxCRWQX2ex34DGgiIj+JSD/vjzJEZKWIfAG4gf+E+C0YY0zMs1pYxhhjjshqYRljjAmomDoDEZEtwPfeZgKwNYzhhIO959gQa+851t4vhP49n6yqf7uMNaYSSEEisvhwp2TRzN5zbIi19xxr7xci5z1bF5YxxpgSsQRijDGmRGI5gTwT7gDCwN5zbIi19xxr7xci5D3H7BiIMcaY4xPLZyDGGGOOgyUQY4wxJRJzCUREuojIGhFZKyJ3hzueUChqZcdoJSINRcQjIl96V768NdwxBZuIVBSRRSKywvueHwp3TKEiInEiskxEYqKoqohs8JZyWi4iYS2tEVNjICISB3wD/AP4CcgFrlLVL8MaWJCJyPnATuAlVW0e7niCzVtss76qLhWRqsAS4OJo/ncWZ3W2yqq6U0TKAQuAW1V1YZhDCzoRGQy0Bqqpao9wxxNsIrIBaK2qYZ88GWtnIG2Ataq6XlX3AZOBQKxPEtFUdT7we7jjCBVV3aiqS72P/wS+Ak4Ib1TBpY6d3mY57y3q/zoUkROB7sBz4Y4lFsVaAjkB+LFA+yei/Isl1nmXRD4L+DzMoQSdtytnOfArMEdVo/49A+OBO4H8MMcRSgrMFpElItI/nIHEWgIxMUREqgBvArep6o5wxxNsqpqnqsnAiUAbEYnq7koR6QH8qqpLwh1LiLVX1bOBrsBAbxd1WMRaAvkZaFigfaJ3m4ky3nGAN4FXVXVauOMJJVXdBniALmEOJdjaAT29YwKTgY4i8kp4Qwo+Vf3Ze/8r8BZO13xYxFoCyQWSROQUESkPXAlMD3NMJsC8A8rPA1+p6mPhjicURKS2iFT3Po7HuVDk67AGFWSqOlRVT1TVRJz/yzmqek2YwwoqEansvTAEEakM/JMwrsgaUwlEVQ8Ag4APcAZWs1V1dXijCr4jrOwYrdoB/8L5i3S599Yt3EEFWX3A412lMxdnDCQmLmuNMXWBBSKyAlgEzFTV98MVTExdxmuMMSZwYuoMxBhjTOBYAjHGGFMilkCMMcaUiCUQY4wxJWIJxBhjTIlYAjHGGFMilkCMMcaUiCUQY2KAiDwhIktFJCXcsZjoYQnEmCjnLXlRB7gBiPr1MkzoWAIxpphE5EERGeJ9/OkR9qsuIjeFLrIjU9VdOKVOPgQmhDcaE00sgRhTAqp63hF+XB2ImAQiIrWASsCfwIEwh2OiiCUQY45ARO4VkW9EZAHQpMD2nd77yiIy07sW+SoRuQJ4BGjkLeI4xrvf294FgFb7FgESkUQR+UpEnvVun+2tpIuI/J+IfOF93ZcLHPca79rny0Xkae8yzUdzH/AosBo4I0C/GmMoG+4AjIlUItIKp0x4Ms7/laU466sX1AX4RVW7e5/jwln9sLl3cSeff6vq794EkSsib3q3JwFXqer1IpINXCoiy3C+9M9T1a0iUtP72qcDVwDtVHW/iDwJ9AZeOsJ7SATOAwYD7XESSJHdb8YcC0sgxhStA/CWqu4GEJHDrR2zEhgrIqOBd1X1YxGpcZj9bhGRXt7HDXESxybgO1Vd7t2+BEgEagBTVXUrgKr61rO/AGiFk4AA4nGWrz2Sh4Hhqqoi8hV2BmICyBKIMcdBVb8RkbOBbsDDIjKPQ84IRCQN6AScq6q7ReRDoKL3x3sL7JqHkxSKIsCLqjq0OLGJSDJwCdBeRCZ6j7myOM81pjhsDMSYos0HLhaReO8qcBceuoOINAB2q+orwBjgbJzB6qoFdnMBf3iTR1Og7VGOmwNc7h38xteFBcwDLhOROr7tInKy9/E8ETnhkNcZDfRU1UTvqn0tsTMQE0B2BmJMEVR1qYhMAVbgdBXlHma3FsAYEckH9gMDVPU3EflERFYB7+GMZ9zo7UJaAyw8ynFXi8hI4CMRyQOWAX1V9UsRuQ+YLSJlvMcbKCI/Ao0BX1cXItIRqKSqcwu87mYRqSIiNQt0ixlTYrYioTGlnIg0xxmkHxzuWExssQRijDGmRGwMxBhjTIlYAjHGGFMilkCMMcaUiCUQY4wxJWIJxBhjTIlYAjHGGFMilkCMMcaUyP8D3BXpoxOHAXoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "vibronic_structure = VibronicStructure1DFD(mol, energy_surface)\n",
    "\n",
    "plt.plot(xdata, ydata, 'kx')\n",
    "x = np.arange(min(xdata) - 0.25, max(xdata) + 0.25, 0.05)\n",
    "plt.plot(x, energy_surface.eval(x), 'r-')\n",
    "plt.xlabel(r'distance, $\\AA$')\n",
    "plt.ylabel('energy, Hartree')\n",
    "dist = max(ydata) - min(ydata)\n",
    "plt.ylim(min(ydata) - 0.1 * dist, max(ydata) + 0.1 * dist)\n",
    "for N in range(15):\n",
    "    on = np.ones(x.shape)\n",
    "    on *= (energy_surface.eval(energy_surface.get_equilibrium_geometry()) +\n",
    "           vibronic_structure.vibrational_energy_level(N))\n",
    "    plt.plot(x, on, 'g:')\n",
    "on = np.ones(x.shape)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a partition function for the calculation of heat capacity\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The partition function for a molecule is the product of contributions from translational, rotational, vibrational, electronic, and nuclear degrees of freedom. Having the vibrational frequencies, now we can obtain the vibrational partition function  $q_{\\rm vibration}$ to compute the whole molecular partition function \n",
    "\\begin{equation}\n",
    "q_{\\rm vibration} = \\prod_{i=1} ^M \\frac{\\exp\\,(-\\Theta_{\\nu_i}/2T)}{1-\\exp\\,(-\\Theta_{\\nu_i}/2T} \\, . \n",
    "\\end{equation} \n",
    "Here $\\Theta_{\\nu_i}= h\\nu_i/k_B$, $T$ is the temperature and $k_B$ is the Boltzmann constant.  \n",
    "\n",
    "The single-point energy calculations and the resulting partition function can be used to calculate the (constant volume or constant pressure) heat capacity of the molecules. The constant volume heat capacity, for example, is given by \n",
    "\n",
    "\\begin{equation}\n",
    "C_v = \\left.\\frac{\\partial U}{\\partial T}\\right|_{N,V}\\, ,\n",
    "\\qquad\n",
    "{\\rm with} \\quad\n",
    "U=k_B T^2 \\left.\\frac{\\partial {\\rm ln} Q}{\\partial T}\\right|_{N,V} .\n",
    "\\end{equation}\n",
    "\n",
    "$U$ is the internal energy, $V$ is the volume and $Q$ is the partition function. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we illustrate the simplest usage of the partition function, namely creating a Thermodynamics object to compute properties like the constant pressure heat capacity defined above. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjiUlEQVR4nO3deZQc5Xnv8e/T07NoFm2j0S6hDSRhFgHD6tiAAYd4tw92IBvBi+zEONg39ybmOPfazonPse91vHEcEmLjJcbYMTYxwTEYCBhjgxIJZLSyCK1IQtJom0XTPd393D+qetQaTUvdM9XTPd2/zzl9uraperqmp55537fqfc3dERERGU6s3AGIiEjlUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbzGPEmY2d1mts/M1ucse6+ZbTCzjJl1jnVMIiIyvHKUJL4NXD9k2XrgPcCTYx6NiIjkFR/rA7r7k2a2YMiyTQBmNtbhiIjIKYx5khgtM1sJrARoaWm5aNmyZWWOSERkfFmzZs0Bd+8oZNtxlyTc/S7gLoDOzk5fvXp1mSMSERlfzGx7odvq7iYREclLSUJERPIqxy2w9wJPA0vNbJeZfcDM3m1mu4DLgZ+Z2cNjHZeIiJysHHc33ZRn1f1jGoiIiJyWqptERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPIa8yRhZneb2T4zW5+zbKqZPWJmL4XvU8Y6LhEROVk5ShLfBq4fsuyTwGPufibwWDgvIiJlNuZJwt2fBA4OWfxO4Dvh9HeAd41lTCIiMrxKaZOY4e57wum9wIxyBiMiIoFKSRKD3N0Bz7fezFaa2WozW71///4xjExEpPZUSpJ4zcxmAYTv+/Jt6O53uXunu3d2dHSMWYAiIrWoUpLEA8DN4fTNwE/LGIuIiITKcQvsvcDTwFIz22VmHwA+D1xnZi8B14bzIiJSZvGxPqC735Rn1TVjGoiIiJxWpVQ3iYhIBVKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8lCRERCQvJQkREclLSUJERPJSkhARkbyUJEREJC8lCRERyUtJQkRE8lKSEBGRvIpOEmb216UIREREKk/8dBuY2b/mzgIrgC+UKqCxsm7XEV472s+1Z88odygiIhWrkJLEUXd/X/h6L/BoqYIxs9vMbL2ZbTCzj5fqOAB/97ONfOhfVvPzdXtKeRgRkXGtkCTxuSHznypFIGZ2DvAh4BLgfOBtZrakFMcaSGf47a7DGPDxH65l/atHSnEYEZFx77RJwt23ApjZtHD+YIliWQ6scvc+d08BvwTeU4oDbdpzlP6BDJ99x+uImXHfml2lOIyIyLhXTMP13SWLIrAeeIOZtZtZM/AWYN7QjcxspZmtNrPV+/fvH9GB1mw/BMC1Z8/g3LmTWLvz8MijFhGpYsUkCStZFIC7byJoEP8F8BCwFkgPs91d7t7p7p0dHR0jOtaa7YeYPamJWZMmcMG8yWzcfZRE6qRDiYjUvGKShJcsiuwB3L/p7he5+xuBQ8CLpTjOs9sPceEZUwBYMW8yyXSGTXu6S3EoEZFxrWJKEgBmNj18n0/QHvH9qI+x+/Axdh/p56IwSZw/bzIAa3ccivpQIiLj3mmfk8hxe8miOO7HZtYODAAfdffDUR9gw+6jwPHkMGtSE9PbGtUuISIyjIKShJnNc/f1eda9zd0fjCIYd39DFPs5lf3dCSBIDgBmxop5k5UkRESGUWh10yNmtmDoQjN7P/DVSCMqsa6eIEm0tzQOLlsxfzLbuvo43JcsV1giIhWp0CTxP4BfmNmZ2QVmdjvwCeDKUgRWKl29SSY2xWmIH//oy2a2AbBlf0+5whIRqUgFVTe5+3+YWQL4uZm9C/ggwZPRb3T3cdXiu78nwbTWxhOWLe5oBWDLvl4uOmNqOcISEalIBd/d5O6PAbcATwCLgDeNtwQBQXXT0CQxd0ozDXUxlSRERIYotOG6m+A5CQMagWuAfWZmgLv7xNKFGK2uniRLpreesKwuZiyc1qIkISIyRKHVTW2lDmSsHOhJcOmik6uUFk9v0QN1IiJD1NTIdKl0hkN9Ayfc2ZS1aForOw72kUxlyhCZiEhlqqkkcTC8xXVa28lJYvH0FtIZZ8fB3rEOS0SkYtVUkjjQHSaJloaT1mXvcHp5n5KEiEhWTSWJrt7wQbrWYaqbsrfBqvFaRGRQ0UnCzP66FIGMha6esCTRenJJorUxzsyJTUoSIiI5Tnt3k5n9a+4ssIJg3Idx50BP/pIEwKKOFl7Zr+omEZGsQm6BPeruH8zOmNmdJYynpA70JGmoizGxafiPvWBaC/+xbs8YRyUiUrkKqW76HBwf4xr4VOnCKa2ungTtrQ0EzwCebGF7C4f7BtTRn4hI6LRJwt23hpN3h/MHSxpRCXX1Jmkfpj0ia+G0FgC2HlCVk4gIVNjIdKV2oCcx7IN0WQuUJERETlBRY1yXWlfPqUsS86c2EzPYpiQhIgLUWEniUF+Sqc35k0RDPMacKRPY2tU3hlGJiFSuYpLEWIxxXTLpjNOXTNPWVH/K7RZOa2XrAT0rISICxY0nMewY1+NFTyIFQGue21+zFrY3s+1AH+7jvnZNRGTUih1PYnARx8eXGBfjSWSTRFvjqT/ygmkt9CRSHOhJ0jFMR4AiIrWkZsaT6OkvsCQR3uG0ratXSUJEat5I+m4638xuDV/nlSKoUujuHwCCPppOZfBZCXXPISJSXJIws9uAe4Dp4eseM/tYKQKLWneBbRJzJk8gHjO2dilJiIgUVN2U4wPApe7eC2BmXwCeBu6IOrCoZaubTtcmEa+LMb+9Wc9KiIhQfHWTAemc+TTj5PmJQu9ugqAPJz11LSJSfEniW8AqM7s/nH8X8M1IIyqRwYbr05QkILjD6ddbDpDJOLHYuMiBIiIlUVRJwt2/BLwfOBi+bnH3r0QVjJl9wsw2mNl6M7vXzJqi2ne2TaKlobAk0T+Q4bXu/qgOLyIyLhVbksDd1wBrog7EzOYAfwGc7e7HwsGObgS+HcX+e/pTtDbGCyoZLMrp6G/WpAlRHF5EZFwq9u6mTjO738yeNbPnzWydmT0fYTxxYIKZxYFmYHdUO+5JDBRU1QTHe4PddkB9OIlIbSu2JHEP8L+AdUAmykDc/VUz+yKwAzgG/MLdfzF0OzNbCawEmD9/fsH770mkCmq0Bpg1sYnGeEx9OIlIzSv27qb97v6Au2919+3ZVxSBmNkU4J3AQmA20GJmfzR0O3e/y9073b2zo6Oj4P13h9VNhYjFjAXtLWxVSUJEalyxJYlPm9k3gMeARHahu/8kgliuBba6+34AM/sJcAXwvQj2TU8iRVuBJQmABdOa2aKnrkWkxhWbJG4BlgH1HK9uciCKJLEDuMzMmgmqm64BVkewXyBouJ45sfCbpRZMa+HxzftJpTPE64ruvUREpCoUmyQudvelpQjE3VeZ2X3As0AKeA64K6r9F1uSWNzRSjKdYeehY4P9OYmI1Jpi/0X+jZmdXZJIAHf/tLsvc/dz3P2P3T1x+p8qTHAL7KkHHMp15vRWAF7ep8ZrEaldxSaJy4C1ZvZCiW6BLYlMxulJFn53E8BiJQkRkaKrm64vSRQl1jeQxv30nfvlmthUz4yJjby0r7uEkYmIVLaikkRUt7uOtUIHHBpqyfRWtqgkISI1rCZu2+lJFDbg0FBnTm/j5X09Gu9aRGpWTSSJ7hGWJBZPb6U3mWbPEXX0JyK1qSaSRHYsiWLaJACWdKjxWkRqW0FJwsy6zezoMK9uMzta6iBHazRtEqAkISK1q6Crpru3lTqQUhoc37rIksS01gYmN9frDicRqVlFjycRdsR3JjDYx4W7PxllUFE7Pr514Q/TAZgZZ81oY/NeJQkRqU3FjifxQeBJ4GHgs+H7Z6IPK1rZNomWxrqif/bsWRN5YW83mYzucBKR2lNsw/VtwMXAdne/GrgAOBx1UFHrSaSYUF83oo76ls9qoy+ZZvtBdRsuIrWn2Ktmv7v3A5hZo7tvBkrS4V+UehKpEZUiAJbPmgjA5j0V3z4vIhK5YpPELjObDPwb8IiZ/RSo+Kew+xIpWopstM46a0YbMYNNShIiUoOK7Zbj3eHkZ8zscWAS8FDkUUWsN5mmuWFkSaKpvo5FHa1s3KPGaxGpPSO7cgLu/ssoAymlY8k0zQ0jq24CWDazjed2HI4uIBGRcaLYu5u+E1Y3ZeenmNndkUcVsd5kalRJYvmsibx6+BhHjg1EGJWISOUrtk3iPHc/nJ1x90MEdzhVtL5EmpYRVjdBcBssqF1CRGpPsUkiFj5MB4CZTWUUVVZjpTeZonmEdzcBnDt3EgDrdh2JKiQRkXGh2Av83wNPm9mPwvn3Ap+LNqTo9SVHV5KY1trI3CkTWLvzcHRBiYiMA8Xe3fRdM1sNvClc9B533xh9WNHqTYyuJAFw/rzJrFXjtYjUmKL/vQ6TQsUnhqx0xkmkMjTXj65W7IJ5k/nZ83vY352go60xouhERCpb1Y8n0Zcceb9Nuc6fNxmA36rKSURqSA0kiTTAiB+myzpn9iTqYqZ2CRGpKaNKEmY2y8wquu6ldxQ9wOaa0FDH0hlt/HbX4QiiEhEZH0ZbkvgXYLOZfTGKYEohqpIEwAXzg8brtLoNF5EaMaok4e7XAouAb0UTTvSOJ4nRlSQALl3UTncixcbdeqhORGrDqNskPLAhimBKoTdsuI4iSVy2cCoAT79yYNT7EhEZD4qqgzGzJuDPgd8BHHgKuDM7xsRomNlS4Ic5ixYB/8fdvzKa/fYlgpLESLsKzzV9YhOLOlp45pWDrHzj4lHvT0SkWO7Bbf09iRS9iRQ9iRQ9/Sl6kym6+1P0JtL0JlJ0Z9f3p+hJ5kyH7bSFKvbK+V2gG7gjnP8DgnaJ9xa5n5O4+wvACgAzqwNeBe4f7X6jLEkAXLaonX9fu5tUOjOike5EpDZkL+Z9yTR9yRTHkmn6kmmODaQHp/uSKY4NZKfTHEumwvdgu95kcMHvTYQJIBlc6FMFtou2NNTR0hintTFOa1OcloY4c6c0F/U5ik0S57j72Tnzj5tZKR6suwbY4u6jHtCoL3t3UwQN1wCXL2rn+6t2sGH30cFnJ0Sksrk7qYzTP5AmkcqQSGWC6YEMiVSa/vB9cHm4TSI7PXDyuv6cZbkX9uzF/thAmmLvcZlQX0dzQx0TGrLvcVob65jW2kxLY5y2xjgt4astvOgPTjcG27Y21tPSWEdLQ5xYzIY9zjf/tPCYir1yPmtml7n7MwBmdimwush9FOJG4N7hVpjZSmAlwPz580+7o95sw/Uob4HNunRRtl2iS0lCxoS7k3HIuJNxxweng3fPHF+X8QK2z10f/izk2SYzgn0Obn98/ydtn/1c4f7TGSeZzjCQzpBKOwPpzAnzwbQzkMqQymRIhtMD6QwDmZzp7HbDTI/mpkQzaIzHaKqvozEeozFeR1N98N4Yj9HWFGfGxEaaG+LBBT682DcNTsdzLvzBfHNDHRPqjyeEpnhd3ot6ORWbJC4CfmNmO8L5ecCLZraOoA37vNEGZGYNwDuA24db7+53AXcBdHZ2nvbXfiyZpi5mNERUNTS9rYllM9t4fPM+PnKl2iXGSiqdOaHo3ZtM05dI0Z9Kkwz/s0umggtLMpXzCucTqeyFwklnnHR48UpnnHR4sUpnPGcZg8vSgxez4S6YORdEH3pBDC6SRW0frs+9qNaauphRX2fUx2LUx2PBdF0sfJ043RCP0TLMuoa4EY+Fy+PBvgYv6vUxmsL3xsH3PAkg3La+zjCrvAv4WCgoSZjZEmAGcP2QVfOAvUCUo/H8HvCsu78Wxc6yAw5F+Qu+Zvl0/vGXr3Ckb4BJzfWR7bdWpNIZ9vck2HOkn9eO9HOgJ8HB3gEO9SU53JfkUF92emCwcS6RyozoWGbQUBejIR6joS5GLGbUmVEXM2IxqDM7cdngOqPOGFzWEI8RC7eNGcG0gVnuvGE564L5nPWx02/PqX5+8Hi563O2jxWzfXZ9Idsc339h+xzmM2Z/luG3icUYvMDXVeB/07Ws0JLEV4Dbh7YRmNlE4Mvu/vYIY7qJPFVNIzHaAYeGc83yGXz98S088eI+3rliTqT7rgbuzr7uBK/s72VbVy/bDvSyvauPPUeOsfdoP/u7E8P+h9zaGGdKSz1TmhuY3NzAgvYW2pqCRrfmhnhQz9p4vO61uSFOU33dYBJojB9PBg3hdDxWu/8BikSh0KvnDHdfN3Shu68zswVRBWNmLcB1wIej2udoBxwazoq5k5nW2sCjm5Qk+pIpNu3pZtOeo4OvF/Z2D7YFQfCf/LypE5g9eQJnzWhj1qQmZkxqYubEJmZOaqKjrZHJExpoiOtuMZFKU2iSmHyKdRMiiAMAd+8F2qPaH4x+wKHhxGLG1Uun8/CGvQykM9TX0K2we44cY/W2Q6zedpDV2w+xac/RwVJBW1Oc5bMmcsNFc1kyvZUF01pY0N7C7MkTVIUgMk4VevVcbWYfcvd/zl1oZh8E1kQfVnR6EykmRPSMRK7rzp7Bj9bs4tcvH+CqpdMj33+l6EumeHpLF0+8sJ8nXtzHzoPHgOBWvQvmT+bWq5dw7tzJLJ/VxpzJE1S1I1JlCk0SHwfuN7M/5HhS6AQagHeXIK7IHBtI097SEPl+r1zawaQJ9dz/3KtVlyT2dffz83V7eXTTa6zaepBkKsOE+jpev6SdW65YSOeCKSyfNbGmSlAitaqgJBHeaXSFmV0NnBMu/pm7/2fJIotIbyLFvKnFPWFYiMZ4HW87bxY/fnYXPYkUrRF0+1FOR/oGeGjDHh747W6e3tJFxmFxRwt/ctkZXLV0OhcvnEJjPPoSmYhUtmLHuH4ceLxEsZRE0CZRmovbey6cwz2rdvDQ+r3ccNHckhyjlDIZ59dbDnDPMzt4bPNrDKSdBe3N3Hr1Et5+/mzOnNFW7hBFpMzG97+/BehNpCIZS2I4F86fwhntzdy3Zue4ShKHepP8aM1Ovr9qB9u6+pja0sCfXL6Ad66YzblzJqldQUQGVX2S6EumRz0qXT5mxk2XzOfzP9/Mxt1HOXv2xJIcJyov7+vhrie38G9rd5NMZbh4wRQ+fu1Z/N65M1WVJCLDquokkUxlSGW8ZCUJgJsuns/XHnuJbzz1Cl9634qSHWc01u48zD8+sYWHN+6loS7G+zrn8keXncGymZWd1ESk/Ko6SfRF3E34cCY11/O+znncs2o7f339MmZMbCrZsYr19JYu7vjPl/jNli4mNsW59eol3HzFAqa1VvSw5CJSQao6SWSf+o36Ybqhbnn9Ar779DbufGILn3nH60p6rEI8t+MQX/zFC/z65S5mTGzkU29Zzk2Xzh/3d2CJyNir6qtGdiyJqLvlGOqM9hZ+/+L5fO+Z7dx8xQIWTmsp6fHy2bj7KF965AUe3bSP9pYG/uaty/mjy86gqV7tDSIyMlWdJAbHkihhdVPWJ647k5+ufZXP/3wT//THnSU/Xq4t+3v48iMv8uDze2hrivM/33wWt7x+YSRDtopIbavqq0hPf1CSaG0sfXfe09ua+LMrF/P3j7zIwxv28ruvm1nyY756+BhfffRF7luzi6b6Oj569WJWvmGxui8XkchUd5LIDl1a4uqmrA9fuZiHNuzlkz9+ngvmTWZ6iRqxD/Qk+PrjL3PPM8HYTzdfsYA/v2oJHW1qkBaRaFV1kuhNZEsSY/MxG+IxvnrjCt76tae49d7n+O77L4m0PeBo/wDfePIVvvnUVo4NpLnhorncdu1ZzJkcWUe8IiInqO4kkRzbJAGwZHob//eG87jtB2u57QfP8fU/uJD4KDvCO3JsgO89s51//tUrHO4b4K3nzuIT153FkumtEUUtIjK8qk4S3f3Z6qax/ZjvXDGHg71JPvvvG3n/d1Zzx40XjKidYO+Rfr71663cs2oHPYkUVy3t4C+vW8q5cyeVIGoRkZNVdZLoTaSIx4zGMox4dsvrFzKhvo7//dP1vPWOX/E3bz2b333djNP2i5RIpXls0z5+tHonv3xxPwBvO282H75yEa+breQgImOr6pNEa1O8bB3W3XjJfM6c0crtP1nHR763hmUz23jXBXPoPGMK89ubmVBfR08ixe7D/Wzcc5RnXuniqZcOcOTYADMnNvGRKxdz48Xzmd8efVfnIiKFqOok0Z1Ilfxp69O56Iyp/MdfvIH71uziB/+9k8//fHPebae3NXLt8hm8/fxZvOHMDg35KSJlV9VJordCBgOK18W48ZL53HjJfPYe6WfD7iPsPnyM/oEMLY1xprc1snz2RGZPalI33SJSUcp/BS2h3kTpugkfqZmTmpg5qXI6ARQROZWqHqS4J5GitUlPH4uIjFT1J4kKK0mIiIwnVZ0keiug4VpEZDyr6iTRk0ipJ1QRkVGo2iTh7vQmUrQ1KUmIiIxU1SaJYwNpMj72XXKIiFSTikoSZjbZzO4zs81mtsnMLh/pvo53E64kISIyUpV2Bf0q8JC732BmDcCI+6PIDjjUpiQhIjJiFXMFNbNJwBuBPwVw9ySQHOn+ehPB0KUqSYiIjFwlVTctBPYD3zKz58zsG2bWMtKdjfWodCIi1aiSkkQcuBC4090vAHqBTw7dyMxWmtlqM1u9f//+vDsb61HpRESqUSUliV3ALndfFc7fR5A0TuDud7l7p7t3dnR05N1Zj5KEiMioVUyScPe9wE4zWxouugbYONL9KUmIiIxepV1BPwbcE97Z9Apwy0h31KtbYEVERq2irqDuvhbojGJfPYkUZtDcoIZrEZGRqpjqpqj1JFK0NpRv6FIRkWpQtUmiV537iYiMWhUnicoblU5EZLyp2iTRrVHpRERGrWqTRK9GpRMRGbWqThIalU5EZHSqNkkc7E0ypbmh3GGIiIxrVZkk3J2DvUnaW5UkRERGoyqTxNFjKVIZp721sdyhiIiMa1WZJA70JgCYppKEiMioVGWS6OoJxiqa2qIkISIyGlWaJIKSRHuLqptEREajKpPEgd6gJKHqJhGR0anKJJEtSUxRdZOIyKhUZZI42Jtk0oR66uuq8uOJiIyZqryKdvXoGQkRkShUZZI40JNgmhqtRURGrSqTRJeethYRiUR1JomehJKEiEgEqi5JpNIZDvUNMFXVTSIio1Z1SeJQ3wCgZyRERKJQdUmiq1dPW4uIRKX6kkTYb5PaJERERq/qksSBHvUAKyISlapLEsd7gFV1k4jIaFVdkli/+whTWxqYPKG+3KGIiIx7VZUk3J1ntnRx+aJ2YjErdzgiIuNeVSWJ7V197D7Sz+WL28sdiohIVYiXO4BcZrYN6AbSQMrdO4v5+d9s6QLgCiUJEZFIVFSSCF3t7gdG8oO/2XKAmRObWDitJeqYRERqUtVUNx3uS/L0li4uX9yOmdojRESiYO5e7hgGmdlW4BDgwD+5+13DbLMSWBnOngOsH7sIK9I0YEQlryqicxDQedA5gMLOwRnu3lHIziotScxx91fNbDrwCPAxd3/yFNuvLrbdotroHOgcZOk86BxA9Oegoqqb3P3V8H0fcD9wSXkjEhGpbRWTJMysxczastPAm1FVkohIWVXS3U0zgPvDRuc48H13f+g0P3NSm0UN0jnQOcjSedA5gIjPQUW1SYiISGWpmOomERGpPEoSIiKS17hMEmZ2vZm9YGYvm9knyx1PqZjZPDN73Mw2mtkGM7stXD7VzB4xs5fC9ynhcjOzr4Xn5Xkzu7C8nyA6ZlZnZs+Z2YPh/EIzWxV+1h+aWUO4vDGcfzlcv6CsgUfIzCab2X1mttnMNpnZ5bX2XTCzT4R/C+vN7F4za6qF74KZ3W1m+8xsfc6yon/3ZnZzuP1LZnZzIcced0nCzOqArwO/B5wN3GRmZ5c3qpJJAX/p7mcDlwEfDT/rJ4HH3P1M4LFwHoJzcmb4WgncOfYhl8xtwKac+S8AX3b3JQQPYH4gXP4B4FC4/MvhdtXiq8BD7r4MOJ/gfNTMd8HM5gB/AXS6+zlAHXAjtfFd+DZw/ZBlRf3uzWwq8GngUoLHCz6dTSyn5O7j6gVcDjycM387cHu54xqjz/5T4DrgBWBWuGwW8EI4/U/ATTnbD243nl/A3PCP4E3Ag4ARPFEaH/qdAB4GLg+n4+F2Vu7PEME5mARsHfpZaum7AMwBdgJTw9/tg8Dv1sp3AVgArB/p7x64iaAnC4bbLt9r3JUkOP5FydoVLqtqYVH5AmAVMMPd94Sr9hLcPgzVe26+AvwVkAnn24HD7p4K53M/5+A5CNcfCbcf7xYC+4FvhdVu3wifJ6qZ74IHD9t+EdgB7CH43a6h9r4LWcX+7kf0nRiPSaLmmFkr8GPg4+5+NHedB/8SVO19zGb2NmCfu68pdyxlFgcuBO509wuAXo5XLwA18V2YAryTIGHOBlo4uQqmJpXydz8ek8SrwLyc+bnhsqpkZvUECeIed/9JuPg1M5sVrp8F7AuXV+O5eT3wjnCskR8QVDl9FZhsZtmHQXM/5+A5CNdPArrGMuAS2QXscvdV4fx9BEmjlr4L1wJb3X2/uw8APyH4ftTadyGr2N/9iL4T4zFJ/DdwZnhHQwNBw9UDZY6pJCx4/PybwCZ3/1LOqgeA7J0JNxO0VWSX/0l4d8NlwJGc4ui45O63u/tcd19A8Lv+T3f/Q+Bx4IZws6HnIHtubgi3H/f/Xbv7XmCnmS0NF10DbKSGvgsE1UyXmVlz+LeRPQc19V3IUezv/mHgzWY2JSyVvTlcdmrlbowZYQPOW4AXgS3Ap8odTwk/5+8QFCGfB9aGr7cQ1Ks+BrwEPApMDbc3gju/tgDrCO4CKfvniPB8XAU8GE4vAv4LeBn4EdAYLm8K518O1y8qd9wRfv4VwOrw+/BvwJRa+y4AnwU2E/Tr9i9AYy18F4B7CdphBghKlR8Yye8eeH94Pl4Gbink2OqWQ0RE8hqP1U0iIjJGlCRERCQvJQkREclLSUJERPJSkhARkbyUJGRcM7N2M1sbvvaa2as58w3lji+XmV1lZleM0bG2mdm0cPoiM9tqZheMxbGlulTS8KUiRXP3LoLnBzCzzwA97v7FcsVjZnE/3o/QUFcBPcBvItpfIT9/HsHT2b/v7s+NdD9Su1SSkKoT/uf8SzNbY2YP53Rd8ISZfdnMVlswHsPFZvaTsG/9vwu3WWDBeA33hNvcZ2bNBez3K2a2GrjNzN4ejl/wnJk9amYzwg4aPwJ8IizlvMHMvm1mN+TE3RO+X2VmvzKzB4CNFoyl8f/M7L/D8QE+XOCpWE7w0N0fu/t/RXJypeYoSUi1MeAO4AZ3vwi4G/hczvqku3cC/0jQjcFHgXOAPzWzbA+hS4F/cPflwFHgz8M+tE613wZ373T3vweeAi7zoCO+HwB/5e7bwmN+2d1XuPuvTvM5LgRuc/ezCJ6uPeLuFwMXAx8ys4UFnIufAre6+1MFbCsyLFU3SbVpJLjoPxJ070MdQXcGWdl+vtYBGzzsz8jMXiHo/OwwsNPdfx1u9z2CgW4eOs1+f5gzPRf4YVjSaCAYB6JY/+Xu2Z97M3BeTqljEsGAMqfb76PAB83sYXdPjyAGESUJqTpGcPG/PM/6RPieyZnOzmf/Hob2VeMF7Lc3Z/oO4Evu/oCZXQV8Js/PpAhL82YWI0gow+3PgI+5++k7YzvRrQSll38ACq2iEjmBqpuk2iSADjO7HIKu1s3sdUXuY37254E/IKg+eqGI/U7ieBfMueMIdwNtOfPbgIvC6XcA9Xn29zDwZ2GVF2Z2lgUDDmFmm0/xOTJh/MvM7G9PsZ1IXkoSUm0yBN1Cf8HMfkvQc26xt52+QDCe+CaCnlbvdPdkEfv9DPAjM1tDMGRm1r8D7842XAP/DFwZ7u9yTiw95PoGQZfYz5rZeoJhJ+PhLa52qg/i7v0ECegdZvbRU39skZOpF1iRHOFdSA+6+znljuV0LBi1b5G7f63csUj1UpuEyDjl7g+WOwapfipJiIhIXmqTEBGRvJQkREQkLyUJERHJS0lCRETyUpIQEZG8/j9fmzcvEoPidAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Q = DiatomicPartitionFunction(mol, energy_surface, vibronic_structure)\n",
    "\n",
    "P = 101350  # Pa\n",
    "temps = np.arange(10, 1050, 5)  # K\n",
    "\n",
    "mol.spins = [1/2,1/2]\n",
    "\n",
    "td = Thermodynamics(Q, pressure = 101350)\n",
    "td.set_pressure(101350)\n",
    "temps = np.arange(10, 1500, 5)\n",
    "ymin = 5\n",
    "ymax = 11\n",
    "\n",
    "plt.plot(temps,\n",
    "         td.constant_pressure_heat_capacity(temps) / const.CAL_TO_J)\n",
    "plt.xlim(0, 1025)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.xlabel('Temperature, K')\n",
    "plt.ylabel('Cp, cal mol$^{-1}$ K$^{-1}$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we demonstrate how to access particular components (the rotational part) of the partition function, which in the H2 case we can further split to para-hydrogen and ortho-hydrogen   components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = Q.get_partition(part=\"rot\", split=\"eq\")\n",
    "para = Q.get_partition(part=\"rot\", split=\"para\")\n",
    "ortho = Q.get_partition(part=\"rot\", split=\"ortho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now plot the constant volume heat capacity (of the rotational part) demonstrating how we can call directly the functions in the 'thermodynamics' module, providing a callable object for the partition function (or in this case its rotational component). Note that in the plot we normalize the plot dividing by the universal gas constant R (Avogadro's number times Boltzman's constant) and we use crossed to compare with experimental data found in literature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# REFERENCE DATA from literature\n",
    "df_brink_T = [80.913535,135.240157,176.633783,219.808499,246.226899]\n",
    "df_brink_Cv = [0.118605,0.469925,0.711510,0.833597,0.895701]\n",
    "\n",
    "df_eucken_T = [25.120525, 30.162485, 36.048121, 41.920364, 56.195875, 62.484934, 72.148692, 73.805910, 73.804236, 92.214423,180.031917,230.300866]\n",
    "df_eucken_Cv  = [0.012287,0.012354,0.008448,0.020478,0.032620,0.048640,0.048768,0.076678,0.078670,0.170548,0.667731,0.847681]\n",
    "\n",
    "df_gia_T = [190.919338,195.951254,202.652107,204.292585,209.322828,225.300754,234.514217,243.747768]\n",
    "df_gia_Cv = [0.711700,0.723719,0.749704,0.797535,0.811546,0.797814,0.833793,0.845868]\n",
    "\n",
    "df_parting_T = [80.101665, 86.358919,185.914204,239.927797]\n",
    "df_parting_Cv = [0.084730,0.138598,0.667809,0.891634]\n",
    "\n",
    "df_ce_T = [80.669344,135.550569,145.464190,165.301153,182.144856,203.372528,237.993108,268.696642,294.095771,308.872014]\n",
    "df_ce_Cv = [0.103048,0.467344,0.541364,0.647315,0.714078,0.798258,0.891147,0.944848,0.966618,0.985486]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABfeElEQVR4nO3deVxUVRvA8d9hR3HHfd8XVFDBNRXNXEuzLLdMWiw1fbVSy7JEszSzTE0rLXczyy0zzTIVM1dQxCUXBCR3BBGRdeC8f5wBUdllmAHO9/O572x37n1mXpuHc+45zxFSSjRN0zTN0liZOwBN0zRNS49OUJqmaZpF0glK0zRNs0g6QWmapmkWSScoTdM0zSLZmDuAnHJ2dpa1atUydxiapmlaFvz8/G5KKcvn9v0FLkHVqlULX19fc4ehaZqmZUEIcfFR3q+7+DRN0zSLpBOUpmmaZpF0gtI0TdMskk5QmqZpmkXSCUrTNE2zSDpBaZqmaRZJJyhN0zTNIukEpWmaplkknaA0TdM0i6QTlKZpFufatWsMGjSIunXr0qpVK3r37s25c+dMcq49e/awf//+dF9bvnw55cuXx83NLXU7ffp0rs6zZcsWZs2aBYC3tzdz5swB4MMPP2Tnzp2AqpRz8+bNHB2rMCtwpY40TSvcpJT079+f4cOH8+OPPwJw/Phxrl+/ToMGDXJ1TIPBgI1N+j93e/bswcnJifbt26f7+sCBA/nqq69ydd60+vbtS9++fR96fvr06Tk6jsFgyPBYhY1OUJqmZWj8+PH4+/vn6THd3Nz48ssvM3x99+7d2NraMnLkyNTnXF1dARg0aBDDhg2jT58+AHh5efHkk08yYMCAh46zfPlyNm7cSHR0NElJSWzatImXX36ZoKAgihUrxuLFiylZsiTffPMN1tbWrF69mgULFtCxY8csP4OUkrFjx/Lnn39SvXp17OzsePnllxkwYEBqvVBnZ2d8fX2ZMGECe/bsYfny5fj6+j6U7B78DLNnz2b79u04Ojryww8/UK9ePby8vHBwcODYsWN06NCB5s2bpx7rwfc7OTkRHR3Nnj17mDp1KqVLl+bEiRM8//zzNGvWjHnz5hEbG8vmzZupW7dulp/VnHSC0jTNopw8eZJWrVql+9rAgQP56aef6NOnDwkJCfz11198/fXXGR7r6NGjBAQEULZsWcaOHUuLFi3YvHkzu3bt4sUXX8Tf35+RI0fi5OTEhAkT0j3GunXr2LdvX+rjAwcOsH37ds6ePcvp06e5fv06TZo04eWXX360D25UqlQpTpw4wcqVKxk/fjxbt24F4NKlS+zfvx9ra2uWL1+erWMdP36cf//9l7Jly1KnTh1effVVDh8+zLx581iwYEGmfyhYAp2gNE3LkKX9gPXq1Ytx48YRHx/P77//TqdOnXB0dMxw/yeeeIKyZcsCsG/fPjZs2ABA165dCQ8PJyoqKstzptfFt3fvXgYPHoy1tTVVqlSha9euj/Cp7jd48ODU2zfffDP1+eeeew5ra+scHcvDw4PKlSsDULduXbp37w5As2bN2L17dx5FbDp6kISmaRbFxcUFPz+/dF9zcHDA09OTHTt2sG7dOgYOHJjpsYoXL26KEDNlY2NDcnIyAHFxcTl+vxAi3fsZfZa050tOTiYhISH1NXt7+9T7VlZWqY+trKwwGAw5ji2/6QSlaZpF6dq1K/Hx8SxevDj1uYCAAP7++29AtWiWLVvG33//Tc+ePbN93I4dO7JmzRpADYxwdnamZMmSlChRgjt37uQoxk6dOrFu3TqSkpK4evXqfa2RWrVqpSbYlBZbTqxbty71tl27dlnun/Z8W7ZsITExMcfntFQ6QWmaZlGEEGzatImdO3dSt25dXFxcmDx5MpUqVQKge/fu+Pj40K1bN+zs7LJ9XG9vb/z8/GjevDnvvvsuK1asAOCpp55i06ZNuLm5pSbBtNatW3ffMPP9+/fTv39/6tevT5MmTXjxxRfvSyRTp05l3LhxuLu757hLDuDWrVs0b96cefPmMXfu3Cz3HzFiBD4+Pri6unLgwAGztBpNRUgpzR1Djri7u0u9oq6maZYks9GERZkQwk9K6Z7b9+sWlKZpmmaR9Cg+TdMKtB07dvDOO+/c91zt2rXZtGlTvsWQ3WHfWs6YrItPCFEdWAlUBCSwWEo574F9PIFfgGDjUxullJlOq9ZdfJqmaQXDo3bxmbIFZQDellIeFUKUAPyEEH9KKR8sZPW3lPJJE8ahaZqmFUAmuwYlpbwqpTxqvH8H+BeoaqrzaZqmaYVLvgySEELUAloAh9J5uZ0Q4rgQYrsQwiWD978mhPAVQviGhYWZMlRN0zTNQpg8QQkhnIANwHgp5YN1RY4CNaWUrsACYHN6x5BSLpZSuksp3cuXL2/SeDVN0zTLYNIEJYSwRSWnNVLKjQ++LqWMklJGG+9vA2yFEM6mjEnTNMtnietBNWnShCVLlpgkBi19JktQQhWR+h74V0r5RQb7VDLuhxCitTGecFPFpGma5UtZD8rT05MLFy7g5+fHzJkzuX79eq6PmVnducwSFKjSSv7+/uzZs4f33nsv23EkJSXlOE7tfqYcxdcBGAacEEL4G597D6gBIKX8BhgAjBJCGIBYYJAsaKUtNK0QGz8e8ng5KNzcILMi6Za6HlSFChWoW7cuFy9exNvbmyNHjhAbG8uAAQOYNm0aoOriDRw4kD///JNJkyZx584dFi9eTEJCAvXq1WPVqlUUK1Ysd19cEWSyBCWl3AeILPb5Cnj0pSo1TSs0LG09qBRBQUEEBQVRr149Pv74Y8qWLUtSUhKPP/44AQEBNG/eHIBy5cpx9OhRAMLDwxkxYgQAU6ZM4fvvv2fs2LG5+VqKJF1JQtO0DFnYclBmWQ8qZcFCe3t7vv32W8qWLcs333zD4sWLMRgMXL16ldOnT6cmqLRLgJw8eZIpU6YQGRlJdHQ0PXr0eJSPX+ToBKVpmkVxcXFh/fr16b724HpQgwYNyvRYeVHZ+8EFC4ODg5kzZw5HjhyhTJkyeHl53bfuU9pzenl5sXnzZlxdXVm+fDl79ux55HiKEl0sVtM0i2Lp60FFRUVRvHhxSpUqxfXr19m+fXuG+965c4fKlSuTmJiYem4t+3SC0jTNoljaelAPcnV1pUWLFjRq1IghQ4bQoUOHDPf96KOPaNOmDR06dKBRo0bZjlVT9HpQmqZpmkno9aA0TdO0QkkPktA0rUCzhPWgNNPQXXyapmmaSeguPk3TNK1Q0glK0zRNs0g6QWmapmkWSScoTdM0zSLpBKVpmsWxlPWgADZv3kzz5s1p3LgxzZo1Y/PmzZnue/r06dTHnp6e6EFduacTlKZpFsWS1oM6fvw4EyZM4JdffuHff/9ly5YtTJgwgYCAgHTP8WCC0h6NHmauaVrGzLAg1K5du/D29mbv3r0PvWaK9aDatm2LtbU15cuXf2g9qGHDhtGlSxdefvnl1Oe+//579uzZw6pVq/D09MTNzY19+/bRv39/Pv/8c0qVKkWpUqXYsGEDr7zyCm3atGH37t1ERkby/fff07FjR+Li4hg1ahS+vr7Y2NjwxRdf0KVLl9x/pxbqUYeZ64m6mqZZFEtaD+rUqVMPPe/u7s7ChQtTHyckJKR2450/f/6hhGkwGDh8+DDbtm1j2rRp7Ny5k4ULFyKE4MSJE5w5c4bu3btz7tw5HBwccvRdFXY6QWmaljELWxDKHOtBZSXt+k/peeaZZwBo1aoVISEhqbGkLFzYqFEjatasyblz51LXlNIUfQ1K0zSL4uLigp+fX7qvPbgeVFbJ4VHXg2rSpMlDsfj5+eHi4pLtc9jb2wNgbW2d6bUw7WE6QWmaZlEsaT2oCRMmMHPmzNSWT0hICJ988glvv/12uvtnd22ptLGcO3eO0NBQGjZsmO3PUlToBKVpmkWxpPWg3Nzc+PTTT3nqqado1KgRTz31FLNnz8bNzS3dcwwaNIjPPvuMFi1acOHChQxjGT16NMnJyTRr1oyBAweyfPny1JaWdo8exadpmqaZhC4Wq2maphVKehSfpmkFml4PqvDSXXyapmmaSeguPk3TNK1Q0glK0zRNs0g6QWmapmkWSScoTdM0zSLpBPWIwsPDeeWVV/jxxx91GRNNyyOWsh7UmTNnaNeuHfb29syZMyfDY7zyyiu4urrSvHlzBgwYQHR0tElizUsPfm4vLy/Wr19vxogephPUI3r//fdZunQpgwcPpm7dusydOzdPClBqWlFlSetBlS1blvnz56db6TytuXPncvz4cQICAqhRowZfffVVrmNNj5SS5OTkPDuewWDIcqFGS6DnQT2CgIAAlixZwpgxY+jevTtz5szhrbfewtvbm3feeYf33nvP3CFq2iMZ//t4/K/55+kx3Sq58WXPLzN8fffu3dja2jJy5MjU51xdXQHTrAf1zTffYG1tzerVqx9aD6pChQpUqFCB3377LdPPVLJkSUAlktjYWIQQD+3j7e3NhQsXCAwM5ObNm0yaNIkRI0YQHR1Nv379uHXrFomJicyYMYN+/foREhJCjx49aNOmDX5+fmzbto1Zs2Zx5MgRYmNjGTBgANOmTXvoPCEhIbz88svcvHmT8uXLs2zZMmrUqIGXlxcODg4cO3aMqlWrsn///vs+N8DevXv54osvuHbtGrNnz2bAgAFIKZk0aRLbt29HCMGUKVOyLNKbV0zWghJCVBdC7BZCnBZCnBJCjEtnHyGEmC+ECBRCBAghWpoqnrwmpWT8+PGULl2aadOm8dRTT+Hj48Phw4fp2LEj77//Pj4+PuYOU9MKnOysBwWkrgeVkqzSc/ToUdavX4+Pjw9Tp06lRYsWBAQE8Mknn/Diiy9Sq1YtRo4cyZtvvom/v/99ySkrvXv35sqVK6mPX3rpJSpVqsSZM2dSl9J4UEBAALt27eLAgQNMnz6dK1eu4ODgwKZNmzh69Ci7d+/m7bffJmV+6vnz5xk9ejSnTp2iZs2afPzxx/j6+hIQEICPj0/qyr4ffvghW7ZsAWDs2LEMHz6cgIAAhg4dyv/+97/U81+6dIn9+/ezcePGdD/31atX2bdvH1u3buXdd98FYOPGjfj7+3P8+HF27tzJxIkTuXr1ara/p0dhyhaUAXhbSnlUCFEC8BNC/CmlTLseci+gvnFrA3xtvLV4mzdvZvfu3Xz11Vep680AeHh48PPPP9OgQQMmTZrEwYMH0/1rStMKgsxaOuZgSetBbdu27b7Hy5YtIykpibFjx7Ju3Tpeeumlh97Tr18/HB0dcXR0pEuXLhw+fJg+ffrw3nvvsXfvXqysrLh8+XJqd2bNmjVp27Zt6vt/+uknFi9ejMFg4OrVq5w+fZrmzZszffr01H0OHDjAxo0bAbUi8KRJk1Jfe+6557C2ts7wMz399NNYWVnRpEmT1Bj27dvH4MGDsba2pmLFinTu3JkjR47Qt2/fXHxrOWOyFpSU8qqU8qjx/h3gX6DqA7v1A1ZK5SBQWghR2VQx5ZX4+HgmTJiAi4sLr7/++kOvOzo6Mn36dA4fPmxxFx01zdJZ0npQOWVtbc2gQYNSE+GDHvxjVQjBmjVrCAsLw8/PD39/fypWrEhcXBxwf/zBwcHMmTOHv/76i4CAAPr06ZO6X3Zld+0qAEuoMpQvgySEELWAFsChB16qCvyX5vElHk5iCCFeE0L4CiF8w8LCTBZndn355ZcEBQUxd+5cbGzSb4S++OKLNG3alPfee4+EhIR8jlDTCi5LWg8qO6SUBAYGpt7fsmULjRo1SnffX375hbi4OMLDw9mzZw8eHh7cvn2bChUqYGtry+7du7l48WK6742KiqJ48eKUKlWK69evs3379nT3a9++PT/++CMAa9asybDbMidrV61bt46kpCTCwsLYu3cvrVu3zvJ9eUJKadINcAL8gGfSeW0r8Fiax38B7pkdr1WrVtKcrl69Kp2cnORTTz2V5b5bt26VgFywYEE+RKZphcfly5flc889J+vUqSObNGkie/fuLc+dOyellDIhIUGWKVNGenl5ZXqMZcuWyTfeeCP1cXh4uOzXr59s1qyZbNOmjTx+/LiUUsqzZ8/KZs2aSVdXV7l37977jnH16lVZtWpVWaJECVmqVClZtWpVefv2bSmllL169ZKXL1+WSUlJsn379rJp06bSxcVFDhkyJHWftKZOnSqHDRsm27ZtK+vVqycXL14spZQyLCxMtm3bVjZt2lR6eXnJRo0ayeDgYBkcHCxdXFzuO8bw4cNl/fr1ZdeuXWX//v3lsmXLpJRSfvDBB/KXX36RUkoZEhIiu3TpIps1aya7du0qL168mPren3/+OfVYD37uB18vXry4lFLK5ORkOWHCBOni4iKbNm0qf/zxx0y/97QAX/kI+cOkxWKFELbGJLRDSvlFOq9/C+yRUq41Pj4LeEopM7wCZ+5isa+++iorV67k1KlT1K9fP9N9pZR07dqVU6dOERgYmDrSR9O0osfb2xsnJ6csh6wXJhZbLFaoztbvgX/TS05GW4AXjaP52gK3M0tO5paYmMgPP/yAl5dXlskJVP/y7NmzCQsLy3SSn6ZpmvYwU47i6wAMA04IIfyNz70H1ACQUn4DbAN6A4FADPDwsBcLcuLECWJjY+natWu23+Ph4cHzzz/P559/zqhRo6hc2eLHgGhagVJQ1oPy9vY2dwgFjl4PKgcWLVrEG2+8QXBwMLVq1cr2+wIDA2ncuDGjRo1i/vz5pgtQ0zTNglhsF19hdODAASpVqkTNmjVz9L569eoxcOBAVq5cSUxMjImi0zRNK1x0gsqBgwcP0rZt21xNvB0xYgS3b9/m559/NkFkmqZphY9OUNkUFhZGYGAg7dq1y9X7O3XqRIMGDViyZEkeR6ZpmlY46QSVTYcOqTnGacuO5IQQgldffZV//vmH06dPZ/0GTdPyzYcffsjOnTtNeo7ly5ffV7svI9lZ9iIkJISmTZtmuc8PP/yQoxgtjU5Q2XTgwAGsra1xd8/19T6GDx+Ora0t3333XR5GpmmWwXuPt7lDyJWkpCSmT59Ot27dTHqe7CaovKITVBFy8OBBXF1dKVasWK6PUaFCBfr168fKlSuJj4/Pw+g0zfym+Ty89ENurV69mtatW+Pm5sbrr79OUlISR44coXnz5sTFxXH37l1cXFw4efIke/bsoVOnTvTp04eGDRsycuTI1LWT/vjjD9q1a0fLli157rnnUhcSrFWrFu+88w4tW7bk559/vq/VUqtWLSZPnoybmxvu7u4cPXqUHj16ULduXb755pvUGD/77DM8PDxo3rw5U6dOBVRSaNy4MSNGjMDFxYXu3bsTGxvL+vXr8fX1ZejQobi5uREbG8v06dPx8PCgadOmvPbaa1nWvvPz88PV1RVXV1cWLlyY+nxISAgdO3akZcuWtGzZMnWNp3fffZe///4bNzc35s6dm+F+Fu1RylCYYzNHqSODwSCdnJzuK5uSWzt27JCAXLt2bR5EpmmWA2/y5DinT5+WTz75pExISJBSSjlq1Ci5YsUKKaWU77//vnz77bfl6NGj5SeffCKllHL37t3S3t5eXrhwQRoMBtmtWzf5888/y7CwMNmxY0cZHR0tpZRy1qxZctq0aVJKKWvWrCk//fTT1HOmLfNTs2ZNuWjRIimllOPHj5fNmjWTUVFR8saNG7JChQpSSvXf8YgRI2RycrJMSkqSffr0kT4+PjI4OFhaW1vLY8eOSSmlfO655+SqVauklFJ27txZHjlyJPWc4eHhqfdfeOEFuWXLlodiSatZs2bSx8dHSilTSw9JKeXdu3dlbGyslFLKc+fOyZTfyN27d8s+ffqkvj+j/UyJRyx1pBcszIZTp04RHR2d6+tPaXXr1o1atWqxZMkSBg0alAfRaZr5eO/xvq/lJKapEa5TO0/F29M7V8f866+/8PPzw8PDA4DY2FgqVKgAqGtFHh4eODg43DensHXr1tSpUweAwYMHs2/fPhwcHDh9+jQdOnQA1PpRaQc5ZVYJPWUpiWbNmhEdHU2JEiUoUaIE9vb2REZG8scff/DHH3/QokULAKKjozl//jw1atSgdu3auLm5AdCqVStCQkLSPcfu3buZPXs2MTExRERE4OLiwlNPPZXuvpGRkURGRtKpUydALaORUiw2MTGRMWPG4O/vj7W1NefOnUv3GNndz5LkKkEJIayAwVLKNXkcj0U6cOAAQK5H8KVlZWXFK6+8wgcffMCFCxeoW7fuIx9T08zF29M7NRGJaQI59dEn/kspGT58ODNnznzotfDwcKKjo0lMTCQuLi51+Yj0lrGQUvLEE0+wdu3adM+T2dITKctOWFlZ3bcEhZWVFQaDASklkydPfmi5nZCQkPv2t7a2JjY29qHjx8XFMXr0aHx9falevTre3t45Xjojxdy5c6lYsSLHjx8nOTkZBweHR9rPkmR6DUoIUVIIMVkI8ZUQoruxZt5YIAh4Pn9CNL+DBw/i7Oyc+hfao3rppZewsrLSgyU0LR2PP/4469ev58aNGwBERESkLkHx+uuv89FHHzF06ND7yhsdPnyY4OBgkpOTWbduHY899hht27bln3/+SV0K4+7du3nWaujRowdLly5NvaZ1+fLl1HgzknZ5i5Rk5OzsTHR0dJaj9kqXLk3p0qXZt28fQOqyIQC3b9+mcuXKWFlZsWrVKpKSkh46X2b7WbKsWlCrgFvAAeBVVC09ATwtpfQ3bWiW4+DBg7Rr1y7PVsatWrUqffr0YdmyZUyfPh1bW9s8Oa6mmdPUzlPz5DhNmjRhxowZdO/eneTkZGxtbVm4cCE+Pj7Y2toyZMgQkpKSaN++Pbt27cLKygoPDw/GjBlDYGAgXbp0oX///lhZWbF8+XIGDx6cOihpxowZNGjQ4JFj7N69O//++29qr4qTkxOrV6/OdLVaLy8vRo4ciaOjIwcOHGDEiBE0bdqUSpUqpXZnZmbZsmW8/PLLCCHo3r176vOjR4/m2WefZeXKlfTs2TO1Zdi8eXOsra1xdXXFy8srw/0sWaa1+IQQJ6SUzYz3rYGrQA0pZe7aonkgv2vxRUREUK5cOT7++GPee++9PDvur7/+St++fdm8eTP9+vXLs+NqWlGzZ88e5syZw9atW80divYAU9fiS0y5I6VMAi6ZMzmZw+HDh4G8uf6UVq9evahQoQKrVq3K0+NqmqYVFlklKFchRJQQ4o4Q4g7QPM3jqPwI0NwOHDiQ2oWQl2xsbBg8eDC//vort27dytNja1pR4unpqVtPhVSmCUpKaS2lLCmlLGHcbNI8LhLLwx48eJBmzZrh5OSU58ceNmwYCQkJ/PTTT3l+bE3TtIIuq1F8fkKIeUKInkIIyx+TmMeSk5M5dOhQnsx/Sk/Lli1p3Lix7ubTNE1LR1ZdfG2ATYAn4COE2CaEGCeEePRhMAXAmTNnuH37dp5ff0ohhGDYsGH8888/BAUFmeQcmqZpBVVWXXwGKeUeKeW7Uso2qKHmd4AZQoijQohF+RKlmaRM0DVVCwpg6NChCCFYvXq1yc6haZpWEOWoWKyU8oqUcqmU8nnAHSjUlSSOHTtGyZIl82TeREZq1KiBp6cnq1atyrJYpKZpWlGSZYISQlQTQkwUQmwRQhwRQuw1tpx6oSbwFlqBgYHUr18/zyboZmTYsGEEBgamrjmlaZqmZT1IYhmwFIgHZgGDgdHATqAnsE8I0cnUQZpLYGAg9erVM/l5nn32WRwcHPRgCU3TtDSyakF9LqXsLqWcL6XcL6UMlFKelFJulFKORQ2eyL8VuPJRYmIiISEh+ZKgSpYsydNPP82PP/5IQkKCyc+naZpWEGQ1SOJkFq8nSCkD8zYkyxAaGkpSUlK+VRsfNmwYERERbNu2LV/Op2maZukyLRYrhDgBpHflXgBSStncJFFZgJQKyPnRggJVfDKl9NHTTz+dL+fUNE2zZFlVM38yX6KwQBcuXADItxaUjY0NQ4YMYdGiRURERFC2bNl8Oa+maZqlyqqL72JmW34FaQ6BgYE4OjpSuXLlfDvn8OHDSUhI4Mcff8y3c2qaplmqbM2DEkK0NQ4xjxZCJAghkgp7sdiUEXymHmKelpubG66urixfvjzfzqlpmmapsjtR9yvUEPPzgCOqosRCUwVlCcy1HLuXlxdHjhzh1KlT+X5uTdM0S5LtShLG0XrWUsokKeUy1DyoQik5OZkLFy7k2wCJtIYMGYKNjQ0rVqzI93NrmqZZkuwmqBghhB3gL4SYLYR4MwfvLXCuXLlCfHy8WVpQFSpUoE+fPqxatQqDwZDv59c0TbMU2U0yw4z7jgHuAtWBZ00VlLnl9xDzB3l5eXHt2jX++OMPs5xf0zTNEmQrQRlH7cVJKaOklNOklG8V1gm6YP4E1bt3b8qVK6cHS2iaVqRldxRfByHEn0KIc0KIoJQti/csFULcEEKkW41CCOEphLgthPA3bh/m5gOYwoULF7C1taV69epmOb+dnR1Dhw7ll19+ISIiwiwxaJqmmVt2u/i+B74AHgM80myZWU7WAyn+llK6Gbfp2YzF5AIDA6lduzbW1tZmi8HLy0vPidI0rUjLboK6LaXcLqW8IaUMT9kye4OUci9QIP/8z68q5plxc3OjefPmejSfpmlFVnYT1G4hxGdCiHZCiJYpWx6cv50Q4rgQYrsQwiWjnYQQrwkhfIUQvmFhYXlw2oxJKc02ByotIQReXl4cPnyY06dPmzUWTdM0c8hugmqDWkH3E+Bz4zbnEc99FKgppXQFFgCbM9pRSrlYSukupXQvX778I542c2FhYdy5c8fsLShQy8Hb2NiwbNkyc4dS5OjFjTXN/LIqFguAlLJLXp9YShmV5v42IcQiIYSzlPJmXp8rJ/K7SGxmKlSowFNPPcXy5cuZMWMG9vb25g6pwIuLg/Pn4dy5e1tQEERGQlQU3LmjtoQEsLMDBwewt1e3ZctClSr3tqpVoX59aNwYKlWCfKyKpWlFQrYSlBCiFDAVSFk91weYLqW8ndsTCyEqAdellFII0RrVmsv0ulZ+MPcQ8weNHDmSTZs2sWHDBoYMGWLucAqchAQ4dAh274Zdu+DAAfVcisqVoV49qF0bSpaEEiXUrZ0dxMerLS4OYmMhPByuXAF/f7h+HZKT7x2ndGlo1AiaNgV3d2jdWt23tc3vT6xphUe2EhRq2feTwPPGx8OAZcAzGb1BCLEWteKusxDiEirB2QJIKb8BBgCjhBAGIBYYJKX5O1YCAwOxsrKiVq1a5g4FgG7dulG3bl2+/vprnaCyyWCAP/6AZctg2zaIiVGtmxYt4H//UwmkQQOVmEqUyP05rl6Fs2fh33/vbRs2wHffqX0cHMDNDdq3h86doWNHKFMmzz6mphV6Ijs5QQjhL6V0y+q5/ODu7i59fX1NdvwXXniBffv2ERISYrJz5NScOXOYOHEiAQEBNGvWzNzhWKx//1VJadUquHYNnJ3h+efhiSdUgsiP5CCl6jI8ckRthw6Br69qiQkBrq7g6Qk9eqiYHB1NH5OmmYsQwk9K6Z7r92czQR0AJkop9xkfdwDmSCnb5fbEuWXqBNW2bVucnJzYuXOnyc6RUzdv3qRatWq88sorLFxYqIvI58rx4+DtDZs3g7U19OkDL70EvXurrjpzi4tTicrHR23796vnHB2hSxcVZ69eUKeOuSPVtLyVXwnKDVgBlEIt9x4BeEkpj+f2xLll6gRVvnx5nnnmGb799luTnSM3hg0bxi+//MKVK1dwcnIydzgW4dQplZjWr4dSpeDNN2HkSKhY0dyRZS42FvbuVd2Pv/0GxnE5NGyoklXv3qo7UI+J0Qq6R01Q2a3F528cDt4caCalbGGO5GRqkZGR3Lx502IGSKQ1atQo7ty5ww8//GDuUMwuLAxefBGaNYMdO+CDDyA4GKZOtfzkBKrl1KMHzJsHgYFqJOGXX0KtWrBokeqSLFcOnnkG1qyB27keiqRpBVumLSghxAtSytVCiLfSe11K+YXJIsuAKVtQfn5+uLu7s3HjRvr372+Sc+SWlBI3NzesrKw4evRovq70aymkhB9/VAMdbt9WLaZJk9SPeWFx964acfjbb/DLL2oghq0tdOsGzz4LffuCiacCalqeMXULqrjxtkQ6W6HrZ7KkOVAPEkIwatQo/P39OXTokLnDyXeXL0O/fjBkiLpWc+wYfPpp4UpOAMWLw5NPwtdfw6VL8M8/KiH/+y+8+qqab9W1K3z1lfpONK0wyzRBSSlTLsTsNC6zkboBf5k+vPyVMgfKEhMUqMoSTk5OfPPNN+YOJV+tXQsuLrBzJ3z+uRpk4JJhYazCw8pKDVGfM0eNDDx6FCZPViMUx46FatWgXbt7r2taYZPdUkcLsvlcgRYYGEilSpUoXrx41jubQYkSJXjhhRdYt24d4eFmn9NscgYDvP22ajU1bQoBAfDWW2qkXlGTMo9rxgw4fVptM2ao4esTJ0Lduur1jz5Sr2laYZBpgjIWh30bKC+EeCvN5g0Uup+JCxcuWOQAibRGjx5NXFwcixcvNncoJnXzJvTsCV98oVoLu3eribWa0rgxvP++alUFBalWVLFi8OGHqnXZuDFMmaK6Qs0//V3TcierFpQd6lqTDfdff4pCVYIoVCxhmY2sNGvWjB49ejBv3jzi4uLMHY5J+Purag/79qmJt/Pn65JBmaldW7U0//lHXZf66itVwmnmTGjZUiX2iRPh4MH7yzNpmqXL6hqUj/F6U9sHrkF9IaU8n08x5ouYmBiuXLlisdef0nrnnXe4fv16oVwrats2dd3FYIC//wYvL3NHVLBUqQJvvKHqDl67BkuWqLJO8+ap61U1aqhBFz4+kJRk7mg1LXPZvQYVY1wPapsQYlfKZtLI8lmQ8SpzQUhQnp6eeHh4MGfOHJIK0a/Mhg3w9NOq6KqfH3hktWazlqny5dXIv+3b4cYNWLlStUyXLFHllqpUUROb//wTEhPNHa2mPSy7CWoNcAaoDUwDQoAjJorJLC5evAhAnQJQb0YIwaRJkwgMDGTTpk3mDidPrF4NAweqH9BduwrGhNuCpHRpGDZMlYMKC4N161SSWr0aundXw9eHD1d/JERHmzlYTTPKboIqJ6X8Hkg0dvu9DHQ1YVz5LjQ0FIAaNWqYOZLs6d+/P/Xr1+fTTz/FAorAP5LFi1VliE6dVBXy0qXNHVHh5uSkiuiuW6eS1ebNqn7hr7/CgAFqblnv3vDNN3qulWZe2U1QKR0AV4UQfYQQLYCyJorJLEJDQ7G1taViAfnT3dramgkTJuDr68vu3bvNHU6uzZ8Pr7+uiqX+9pv68dTyj6OjmgC9cqXqBtyzR13DOnsWRo1Sc608PNTw9ePH9YhALX9lt1jsk8DfQHXU/KeSgLeU8lfThvcwU5U6Gjp0KAcPHkytJlEQxMXFUatWLdzc3Pj999/NHU6OrVypupX691cljCyh8rimSKmqV2zZokouHTqknqtRQw3/79VLVbQoWdLckWqWLF+KxQK3pJS3pZQnpZRdpJStUBXNC43Q0NAC072XwsHBgfHjx7Njxw78/f3NHU6ObN8OL78Mjz+uKkXo5GRZhIAmTeDdd9UqxFeuqIUYW7aEH35Qf1SUK6euY336qW5daaahK0kYhYaGUr16dXOHkWMjR46kRIkSzJ4929yhZNvhw+paR7NmsHGjXlaiIKhUCV55BTZtgvBwNXH67bchMlIlMTc3qFpV/dHx009qorWmPSpdSQIwGAxcvny5wLWgAEqXLs3IkSNZt24dpwtAjZvz59UF+YoVVStKdxEVPHZ2quU0a5aaVH35Mixdqtaw2rRJjcYsX14lrTffVIMv9JIhWm7oShLA1atXSUpKKpAJCmDSpEk4OTnx3nvvmTuUTF27ptZBArWOU6VK5o1HyxtVqqgVjFNGBe7fr+oEliunqrL37Qtly0KbNqq19ccfalkRTcuKTWYvSil9AB8hxHIp5UUhhJPx+UI1U6KgDTF/kLOzM++88w7vv/8+//zzDx06dDB3SA+Jj1fXLa5fVyPF6tc3d0SaKdjYqIoV7dqpWoFxcarE0q5davv8c3XNytZWzXnr0EFt7dtDhQrmjl6zNNm9BlVCCHEMOAWcEkL4CSGamjCufFXQExTAuHHjqFy5Mu+8845FzosaP179UK1YoStEFCUODqo7cPp0VVvx1i34/XfV9QdqmkH//qrLt359Vdpq8WI4dUrXDdSyaEGlsRh4S0q5G0AI4Wl8rr1pwspfKQmqIA6SSFG8eHG8vb15/fXX+fXXX+nbt6+5Q0q1dKma9PnOO2pwhFZ0OTmpbt6Urt64OFXW6p9/VNfgb7+pP2IAypRR3YIeHqq15e6uuhO1oiO786COSylds3ouP5hiHtQbb7zB2rVriYgo2CPnDQYDLi4u2NjYcPz4cWxssvv3h+kcOaIunnfsqAZFWEBImgWTUg2k2b9fJa2DB9X6VimtqUqV7iWrlK2AzK0vkh51HlR2fy6ChBAfAKuMj18ACs0angVxDlR6bGxsmDlzJs8++ywrV67k5ZdfNms8N27As8+qH5W1a3Vy0rImhKq+3qDBvUr2d++qeVa+vmrz81MtrZS/ratWBVdXNW0hZWvUSM+tKwyy24IqgyoS+xggUVUlpkkpb5k2vIeZogXl6upKzZo12bJlS54e1xyklLRr145Lly5x/vx5HB0dzRKHwaCKkB44oP4SbtnSLGFohVR0tFqMMSVhnTihKl+kVGW3sVFJKm3SatIEatYsmisym4vJW1BCCGtgo5SyS25PYulCQ0Pp2LGjucPIE0IIPv30Uzw9Pfnyyy+ZPHmyWeL4+GM1mXP5cp2ctLzn5HSv6zhFYqKqIXjixL3tn39U6z2Fvb0ajNGw4b2tUSN1W6pU/n8OLXNZJigpZZIQIlkIUUpKWeim20VFRREZGVkouvhSdO7cmX79+vHRRx8xcODAfF9CZP9+NWpr2DBVa0/T8oOtLTRtqrbBg+89f/s2nDypWlhnz95LYps3379oY8WKqmuxTh21SnHa20qVwCq7Y561PJPdqwLRwAkhxJ9A6hQ7KeX/TBJVPvrvv/+Agj3EPD0LFiygSZMmjBw5kh07diCEyJfzRkXBCy+orpSvvsqXU2papkqVujffKq2EBAgKupe0zp5VAzT++ktVx0h79cPeXiWrlK1GDVXpvXp1tVWpokt2mUJ2E9RG41boFIY5UOmpXr06M2fOZOzYsaxevZphw4bly3nHjIHQUNi7V5cx0iybnZ3q3mvU6OHX4uLg4kUIDlZJLO3t/v3pl26qWFElq5TEVaWKeq5iRdUCq1hRTUa2tTX9ZyssspWgpJQrTB2IuRTWBAUwatQo1qxZw5tvvknPnj0pX768Sc+3di2sWgXe3qoygKYVVA4O965RpefOHbh0Cf777+Hb8+dV1YyoqPTfW67cw4kr7WNnZ1UaqmxZ1foryl2LmSYoIcSvqAm5v0spEx94rQ7gBYRIKZeaLEITCw0NxdramsqVK5s7lDxnbW3NkiVLaNmyJW+99RarVq3K+k25dPGiWuAupcSNphVmJUpA48Zqy0hMjCrtdf26qkOZcj/t4yNH1P3oDIrHCaEmLKckrMy2UqVUr0XJkiq+kiULfmstqxbUCOAt4EshRAQQBjgAtYFA4Csp5S/pvVEIsRR4ErghpXyoLJJQF0XmAb2BGMBLSnk0tx8kt0JDQ6lWrRrWhXTsadOmTXnnnXeYMWMGL7zwAj1SpvDnoeRktWR7cjKsXq3nO2kaQLFi965ZZSVtMrt5EyIiHt5u3VK3gYH3Hmc1S8jB4f6Eld59JycVa042W1uVPE0tW/OgAIQQtYDKQCxwTkoZk8X+nVCDK1ZmkKB6A2NRCaoNME9K2SarOPJ6HpSnpyfJycns3bs3z45paeLi4nBzcyM+Pp6TJ09SvHjxPD3+woXq2tPSpaqqtaZpppecrK6FRUSoNbpu31Zdj1FR928PPvfg4/j4nJ/b2vrhpOXoqAaK2Nur63v29vDLL/lTSQIpZQgQkoP99xqTWkb6oZKXBA4KIUoLISpLKa9m9xx5ITQ0lPaF/IKJg4MDixcvxtPTkzfeeINly5bl2ai+ixfVEgo9etyb+a9pBZH3Hm+8Pb3NHUa2WVmp7r8yZaBu3dwfx2CA2FjVistqy2y/u3fVyMj4eJUEc5P4HmTOzpiqwH9pHl8yPpdvCSopKYlLly4VygESD+rUqRMffPAB06dPp2PHjrzyyiuPfEwp4bXX1P1vv82fJr+mZSW3iWaaz7QClaAyE58QT1RMFHdi7hAdG82dmDvcjbvL3bi7xCbEEpsQS1xi3ENbvCH+3pYUT4IhgYQktRmSDSQmJZKUnESiTMRgbyDJNglDSQNJyUkkSbUly2SSZRKSZAh4tM9RIK4WCCFeA16DvB1td/36dRITE4tEggL48MMP2b9/P2PGjKFVq1a4ubk90vFWrFCLzy1cqOY9aZolsMREk5ycTGR0JNcjr3Mz6iZhUWFE3Ikg4m6ESiSxd7gTf4fouGjuxt8lJiGGmMQYYg2xxCXFkZCcQIJMwICBRBIxYCBJJJFkpTZpLZFWasOa7C+klEMC4+EFWAuwEeqxjQAbeW9zSAbrPFguJVsJSgjxDPCblDIPGm2pLgNp17eoZnzuIVLKxajRhLi7u+fZYkeFeYh5eqytrVmzZg0tWrRgwIAB+Pn5USqX9V2uXlVr+nTqBCNH5nGgWqFnKd1p3nu8meYzLfWxmKa6AaZ2noq3pzcJCQlcj7jOpbBLXIm4wrXIa1yPus7N6JvcvHuTyNhIbsff5k7iHWIMMcQmxRIv40kgAYOVgSSbJJJtktXa5DlJGlaAjSpdJqwEVklWWCdbY5tshX2yNcWTrXFMtqVYshXFkqxwShIUNwiKJUGxJEmxxGSKJyTjmJiMY0KS2uINOMQl4mCQ2CeBXTY222SBnZ0jtnYOWDsUw8qxGMLBQV1wStkyeSx4tFW+s1ssdhnQFdgLrEMNOzdk4321gK0ZDJLoA4zh3iCJ+VLK1lkdMy8HSfz0008MHDiQEydO0LRpoVl/MUv79u3D09OTfv36sX79+hxfj5ISnnlGLTwXEKBXx9VyTkwTyKl5t7Dmg4kmRUqiSU5O5mbETS5cuUDwjWBCb4Zy+dZlrkVdIywmjIi4CE4UO0G5W+WIJZYEqwSSbJKQ9jLrP+MTwSrRCmuDNTbSBjtph4OVA/bCHkdrR4pZF6O4bXGc7JxwsnOipENJSjmUpJx1MZyTbShrgLLxyZSJNVAmJoGSUbGUvHUH25sREBamhvXduqVGQWT1e12yJJQurbZSpdQwvZSheult6b1WvPi9UQ+POE49X5bbkFK+JISwBXoBg4GFQog/pZSvZhLYWsATcBZCXAKmArbG430DbEMlp0DUMPN8H/9V1FpQKR577DFmzZrFxIkTmTdvHuPHj8/R+9evV3XMZs/WyUkzn+TkZK7duMbp0NM433JmfKXxhN4KZWP8RupF1SPKEMVnGz9jxi8zSLJPgowK+1uBsFZ/pMWJOIpRjPJW5SlhU4LSNqUp41iGcsXK4ezkTPkS5alYuiKVy1SmarmqVHOuRqniaXohEhPVWPErV+7f/ruiuh1uBKmkExamylWkx8EBype/tzVooCY6pSSejLZSpQpdqfZsDzMHMCapnqhk0klK6WyqwDKSly2o//3vf6xcuZLIyMg8OV5BIqWkf//+bN26lY0bN2Z7Bd6oKFUapnJlOHRIz3nSsi+rVk5aMTExBAQG4B/iz9nLZwkKC+JS1CXCYsO4nXybGOsYDI4GKI66MPIAqxgr7A32FKMYJW1KUta+LM7FnKnoVJHKpSpTvVx1ajjXoG7lutSpVAcHW4esux2lVIucXbyottDQe/f/+08lohs3Hm7lWFurEhGVK6taRxUq3J+AHtyKFy80I44etQWV3S6+XsBAVItoD/AT8Ed2uvnyWl4mqKeffpoLFy5w4sSJPDleQXPnzh0ef/xxAgIC2LFjB507d87yPePHw/z5cPiwWs1UK/hSfpgz+oHO6+tFBoMB249t+bndz5wIPcHZ62cJiQzhauxVbiXfIsY2hiSnJHXdJi0JNvE2FEsqRimrUpSzK0el4pWoXro6tZ1r06BKA1xquPDD2R+Y3nV67oKLiVFF986fVzNiAwNVAb6UhPRgq6dECTVCqEYNtXJilSoPb+XLF7qWTXblV4Jai7r2tD2PB0rkWF4mqJYtW1K5cmV+++23PDleQXTz5k06derE5cuX2bNnDy1atMhw32PHVFIaOVKN3NMKh5TrQRldF8rN9aLIyEjOBp7l8NnD+If6czbsLP/d/Y+byTeJcYiBdPperOOtcTI4Uc66HJWLVaZW6Vo0rNQQlxouNK/dnJrlamJrnQe1exITVeL59997iSjl9vID47TKlVPrbdSsef9Wo4a6LV260LR2TMGk16CEEPWAilLKwQ883wG4JqW8kNsTW4LQ0FDatMmyeEWh5uzszB9//EGHDh3o0aMH+/bto0GDBg/tl5ysau05O6vFCDXt7t27nD9/Ht/Tvhy8cJCT104ScjeEcBGOoaQBSnNv5FopsCpuRUlDSarbVcdga6BP9T40rdaUlnVb0qhSI4rb5W2FExITVeI5fRpOnbp3e+7cvaV3QXW51asHjz+uLqrWq3dvK106b2PSciSrKwhfAuktyRplfO2pPI4n39y9e5fw8PAiN0AiPdWqVeOPP/7gscceo3v37uzbt49q1ardt89336lrTqtW6f9mC4OMhlen3Hau2Rmfiz4PvV4vph5cgSuJV4hxjIHyqOtAABXAKsmKssllqeZYjXpl69GsajPa1G9D82rNqeRUyXTrkkVEqCb+0aPqNiDg/kQkhCqK5+ICTz6pbhs3VgMQ9LowFivTLj4hxBEppUcGr52QUjYzWWQZyKsuvjNnztC4cWNWr17N0KFD8yCygs/Pz48uXbrg7OzM9u3baWhca+DGDTUwwtVVLSOgezQKnsyuI6Xt4tvdeTfHTx5n8X+LMQQbCI0PJa5lnFqmNE0Dxy7Jjko2lahfqj5u1dzo0LADLaq1oEapGlgJE64PIaUajJA2GR07pq4RpaheHdzcVBJycYEmTdQ/4GLFTBeXli5TDzMvnclrGQ3aLBCK6hDzzLRq1Yq//vqLPn360L59e3799Vfat2/PpElqOYBFi3RysjTZHcCQtrpCVFQU4zePp01sG46dPAbOUOKJEvAYdFnXRbWKigEuYC3Vxf2edXrSoW4H2tdtj0sFFyoUr5A/qzTHxoKvLxw4oLaDB9X6FKD+Mdavr9Z4GT0aWrZUick53wcXayaSVYLyFUKMkFIuSfukEOJVwM90YZmeTlDp8/Dw4MCBA/Tq1YvHH3+cDz7YzooVnkyenPnaN5p5ZFXW59q1axw7dgyAZ597lsMhh7kkL8FTsOzaMqig9ot+LBo77GjWsBmta7bm64CvOTnqJA2dGzJj74z8qfogpWoJpSSjAwfA319VMwV1TeiJJ8DDQyWj5s3VKDqt0MoqQY0HNgkhhnIvIbmjBoD2N2FcJhcaGoqVlRVVqlQxdygWp27duuzfv58nn+zL+++XoHTpO7z/vv4hyE85HdotpSQoKIhjx46lbj4OPsS43lsVZ2PTjZCmYErHlh3pVLcTrau2xqOKB9/6fcs0n2n4Baj/1Jt+rXae2nlqnnymdIKGCxdg92617dmjJrOC6o7z8IAJE1QLqW1bNZhBK1KyO8y8C/f+aZ+SUu4yaVSZyKtrUF5eXvz111/8999/We9cRH39dTyjR9sDQxg2zIavvvqKkvqCcp7z3uOtbtMkpMyGdk/dPZXpe9OZ53MYuA60BlsnWxKLqwEC9lb2xCfH07ZqWw5ePvjw8dKZKJvXpYhShYbeS0i7dqkJrqAmsnbpAh06qITUvLmeBV4I5Fepo93A7tyexBKFhobq7r1MREXBtGn2tGsneeKJBsyY8RH79u1j9erVhX79rPyWMpouvRaTlJLLly9z5MgRDh8+zJEjR/D19VXjaKsAI6DMjTLEOMcQ3/reFMWnGj/FY9Uf47Eaj+FWyQ27GXYcePVA6usmS0APun0b/vxTlb3ftUu1mEBdJ/L0VIuJde0KDRvqC5zaQ4rsnyihoaG461IIGfr4Y1VSbOtWgbu7Nz16dOeFF16gY8eOTJkyhQ8++AAb/Rdunspo6Dd71GZd0ZrKHSpTZlQZEhwTiJWxAJRvXD41GT1W4zEafNWADc9vuO/YOe2my3W3npRqvtG2bWrbt09dQypVCjp3hrFjVUupaVO14p6mZSJHtfgsQV508SUlJeHg4MCECROYOXNmHkVWeJw/r0bnDh0Ky5bdez4qKoqxY8eycuVKPDw8mDdvHu3atTNfoAVYRnXpAJWQPKHuL3Wp2LYihuoGgkUwYfFhANQqXYtutbvxeJ3H8b3iy5zuc3JU5y7l/Hk28CEmRrWOUpJSypDv5s2hTx/o3VtdQ9J/0BQ5+VLqyJLkRYK6ePEitWrVYvHixYwYMSKPIis8+vVTvzfnzqn6lg/66aefGDduHNeuXeP5559n1qxZ1K5dO/8DLUASEhI4ceIER44cSd1OnTpFcqdkVeESeDrgaSq4VyCyQiQ/Xfop9b3lHMvxeJ3Hebz243Sr0406Zepkeq586b6LiIBfflGl7f/6S63vXbw4dOumklKvXvDAZG+t6MmXa1CFTXBwMAB16mT+H3pR9OefsGULzJqVfnICeP755+nduzefffYZn332GZs3b2bcuHFMnjyZMmXK5G/AFigpKYkzZ87cl4yOHz9OQkICAOXKlcPDw4Onn36aj6w/Sn3fbo/d3I6/je0VW2qXrs1oj9F0q9ON5hWbm3bya3bduKHWWdmwQf0FYzBArVqqOGOfPmr1Snt7c0epFSJFsgW1dOlSXnnlFS5cuKCTVBoGg6oWERenLiNk57fm8uXLTJkyhRUrVuDo6Mjw4cMZO3YsjYvIpCkpJRcuXMDX1zc1GR09epS7d+8CUKJECVq1akVC+wTGuY7D3d2dKMcofjv/G1vPb+XgJTWqrrhtcQa6DKRPgz50q9ONkva5Hy2Zp913V6/Cpk2qpeTjo4oy1qsHAwaorWVLPbhBy5BuQeVCUFAQ1tbWVK9ePeudi5BvvlGJadOm7P8hXLVqVZYtW8Zbb73F3LlzWbp0KV9//TXdu3dn7Nix9OzZs0APpkj7Y5+YmMjZs2f5cNeH1LtUD39/f3x9fbl16xYA9vb2tGjRApc3XBjTdAzu7u40bNgQIQRW063oUKYD7/76LsGRwQ+d527iXaqXqs4zjZ959JgfNTndugU//QRr1qhBDlKqUkHvvaeSUvPmOilp+aJItqCGDh3KgQMHCAoKyqOoCr7wcFU1pmVL1c2X29+fsLAwlixZwsKFC7ly5Qply5alb9++PPPMMzzxxBM4ODjkbeAmcuvWLY4fP04Xny68dPEljh8/zsmTJ1U3nTfYfWKHi4sL7u7ueHh44O7uTtOmTbG1tUVMEyR9mMSB/w6w/vR6Nvy7gf+i/sPGyoZudbrxbONnebLBk6p4an4N985KQgL8/ruqBrxli3rcuDEMGqSSUpMm5o5QK4B0CyoXgoKC9EX9B3h7qykrX375aH8cly9fnvfee4+JEyfy22+/sXHjRjZt2sTy5ctxcnKie/fudOjQgfbt29OiRQvszXjNIiEhgeDgYM6dO8e5c+c4e/Zs6v2rKRUNvOG3337D1dWV//3vf7i5ufFC4AtER0dja3v/2kRJyUn4hKgK4NW+qMbV6Kv3vW5INvB74O+0qdqGSk6V8uMjZk5KVedu5Ur48Ue4eVMtrjdqFAwbprvvNLMrki2oihUr0rdvX5YsWZL1zkXAyZOqxubrr5tmIcKEhAT27NnDhg0b+OOPPwgJCQFUl1irVq1wdXWlfv36qVvt2rWxs3twOdWci4mJ4caNG1y6dCl1+++//7hw4QLnzp0jKCiIpKSk1P2dnZ1p0KABd1vf5Xjp49k6R+eanZnfaz4jt47kwKUDD73+bod3mdltZoYtpbxerTZbLl+GFStUa+nMGdWf268fvPgidO8OtnmwKKCmoYeZ59jdu3dxcnLik08+YfLk9Ja6KlqkVL9Jfn5q/lO5cqY/57Vr1zhw4AD79+9n//79nD59msjIyNTXhRCUKVOGcuXKpW6lSpXCxsbmvi0xMZG4uDhiY2OJjY0lJiaGiIgIwsPDCQ8PJzY29qFzOzk5UbduXRo0aJC6NWzYkPr161O2bNmH9k+bWNIuS3HlrSv8cOIHJvw5AQAbKxt61evF0GZDGbRh0EPJyOxdeUlJqprDt9/Cr7+qwQ4dO6qkNGCAXuRLMwndxZdDeoj5/X79FXbuhPnz8yc5AVSqVIn+/fvTv7+qNyylJDw8nPPnz3P+/HmCgoIICwtLTTZXr17l7NmzGAyG1C0xMRFbW1scHBxwdHRM3WrUqEGLFi1wdnamXLlylC9fnurVq1OtWjWqVatGyZIlH6nVsjpgNQDV5lYjWSYDsKDXAga6DKR88fIADNow6KH3mazgalauXYOlS2HxYjWBtkIFmDQJRoxQS5lrmgUrci2oLVu20K9fPw4dOkTr1q3zMLKCJz5eVYyws4Pjx4tOz05OWjOdl3Vmb+jebO2bUrXBLN12aSUnq8mz336rJtMaDGo589dfV115edB9qmnZoVtQOZQyck+3oGDePFW78/ffi05yyq7IuEhWB6wmIi4CgBJ2JbiTcIe9XnvpUKND6sTZ9JKd2ZJTZCR8952aL3DhgmoSjx8Pr72mhmhqWgFjAdPT81dwcDAlSpSgXH71Z1moa9dgxgx46ino0cPc0Zie9x5vxDSRWoA15X7KUheguhoPXjrIy7+8TJXPqzB2+1jsre1Z/ORirrx9BYCONTtaRlWHtM6fhzFjVGmhiROhShX44Qc1GOKzz3Ry0gqsItmCql27dv4sV23B3n9fVYz4/HNzR5I/vD3vdbs92Oq5HXebNSfW8K3ftwRcD8DJzolhzYfxWqvXaFWlVep+6V1HMtu1JSnVmkpffglbt6om8JAhqsXk6mqemDQtjxXJBNWgQQNzh2FWvr6qSvnbbxftP67Ph59n3qF5LPdfzt3Eu7So1IJv+nzDkGZDKGH/8ArC6VYFz+/uvPh4WLtWJabjx9W8pQ8+UHOXKlnA3CpNy0NFKkFJKQkODqZnz57mDsVspFR/ZJcvD1OmmDua/CelZLjrcPqu7cvWc1uxtbZlSLMhjHYfjXsVd8ttWYeHw1dfwaJFqmhr06bw/feq1VRAqnNoWk4VqQR1/fp1YmNji3QViXXr4J9/1LX0UqXMHU3+SUhK4MeTPzL34Fz8r/njXMyZDzp9wCiPUZZR1SEjly6pftjFi9W6S717w5tvqlF5lppMNS2PFKkEVdRH8MXEqGvoLVqAl5e5o8kfN2Nu8q3vtyw8spCr0VdpUr4JS55awtBmQ3G0dTR3eBk7dw4+/VRVe0hOVqtHvvOOromnFSk6QRUhs2erP8h/+AGsrc0dTc7lZH5R6O1QZv8zm6XHlhJriKVH3R4sf3o5T9R5wnK78QCOHYOZM9XyFvb2aoj4hAlq3SVNK2KKVIJKqSJRqwj+xx4aqhLUwIGqwk1BNM1nWpYJKjAikFn7ZrHi+AoEghddX+Stdm/RpLwFtzykhL17VWLasQNKloR334Vx46BiRXNHp2lmU6QSVFBQEFWqVCkwSz7kpQmqZByzZ5s3DlM5HXaamftm8sOJH7CztmOU+ygmtp9I9VIWvOZXylBxb2/4+29VhuiTT2D06KJ1gVDTMmDSGYdCiJ5CiLNCiEAhxLvpvO4lhAgTQvgbt1dNGU9QUFCR7N7btQt+/hkmT4YaNcwdTc5kNcHW/5o/z/38HE0XNWXTv5t4q+1bBI8LZn6v+ZabnKRU/6d07qwGOwQFwYIFEBKi/k/SyUnTABPW4hNCWAPngCeAS8ARYLCU8nSafbwAdynlmOwe91Fq8dWoUYMuXbqwYsWKXL2/IEpMVIMiYmLUarkFufGYdoLt4cuH+WjvR2w9t5WS9iX5X+v/Ma7tOJyLOZs5yiyktJj27lUVHyZPhldfLdj/x2haBiy5Fl9rIFBKGQQghPgR6AeczvRdJhIfH8+lS5eK3BDzRYvg1CnYvLlw/AaevHGSKbum8MvZXyjrWJaPunzEmNZjKO1Q2tyhZW7PHpWYfHxUYlqwQCcmTcuCKbv4qgL/pXl8yfjcg54VQgQIIdYLIdLtkxFCvCaE8BVC+IaFheUqmIsXLyKlLFJdfDduwNSpqtZe377mjiZjnss9s9wn6FYQzSo0o/nXzdkdspuPunxEyLgQpnSaYtnJae9e8PSELl3U0PH581Uh1zFjdHLStCyYu+rlr0AtKWVz4E8g3b43KeViKaW7lNK9fPnyuTpRURxiPnky3L376Mu4m5rPRZ8MX7ty5wqjfxtNw68aEhgRyMT2EwkeF8yUTlPSLUdkMfz81F8GnTvfS0xBQTB2rE5MmpZNpuziuwykbRFVMz6XSkoZnubhd4DJxpilDDEvKl18hw+rdeomTIBGjcwdTc5FxEbw6b5PWXB4AYnJiYxoOYIpnaZQpUQVc4eWuTNnVG289evVchdz5qhReY4WPClY0yyUKRPUEaC+EKI2KjENAoak3UEIUVlKedX4sC/wr6mCCQoKwt7ensqVK5vqFBYjKUn1IFWqpH4rLZHncs/7Wk4po/Qeq/4YPev1ZPb+2dyJv8PQ5kOZ5jmNOmUsvOUbGgrTpsHy5VCsGHz4oarGW7KkuSPTtALLZAlKSmkQQowBdgDWwFIp5SkhxHTAV0q5BfifEKIvYAAiAC9TxZOyzIaVlbl7NU3v22/hyBFYs8Zyfx/3eO1JvS+mCZI+TGLV8VW8t+s9puyewtONnuajLh/RtEJT8wWZHTduqLlLX3+t+lHHjVN9q7nsitY07R6TTtSVUm4Dtj3w3Idp7k8GJpsyhhTBwcFFonvvyhX1+/jEEzB4sLmjyT6PJR4cvXqU1lVb89OAn+hQo4O5Q8rc7duqiOvcuWoM/0svqREp1S107pWmFUBFopKElJILFy7Qvn17c4dicuPHqyWDFi2y7IERoNZjmrRzEqCKuv7wzA8MbDrQ8lasTSs2Vi17MWsWRETAc8/BRx9Bw4bmjkzTCp0ikaBu3bpFVFRUoR/Bt22bqhgxYwbUq2fuaDIWERvBdJ/pLDyyEAcbBz7p+gnj24637OriiYlq1Mn06aqZ2rMnfPwxtGxp7sg0rdAqEgkqZQRfYU5QMTHwxhvQuLFaUsMSJSQlsOjIIqb7TOd2/G1ebfEq07tMp6KTBRdETU5Wi2h9+CEEBkL79mpF206dzB2ZphV6RSJBpcyBKszXoKZPV6XcfHzAzs7c0Txs+/ntjN8xnnPh5+hetzufd//csgdASAm//Qbvvw8BAdC8OWzdqhYMtPS+U00rJHSCKgROnFDX6195xfL+sA+MCOTNHW+y9dxW6petz29DfqN3/d7mDitze/fCe++ppYfr1lULaA0cCEVgBKimWZIikaBOnTpFpUqVKGmpY64fQWKiGkBWpoxagNVSRCdE88nfn/D5gc+xs7ZjdrfZjGs7DjtrC2zepTh6VCWmHTtUvbxvvoGXXwZbW3NHpmlFUpFIUIcOHaJ169bmDsMkPvlEVdXZsEEVLjA3KSVrT65l4p8TuXLnCi+6vsisx2dRuYQFT5BOW/2hbFn47DN1QU9Xf9A0syr0CerWrVucO3eO4cOHmzuUPOfnp0bsvfACPPOMuaOBY1ePMXb7WP757x9aVW7F+ufW0656O3OHlbG01R8cHVWSevttvR6TplmIQp+gDh8+DEDbtm3NHEneiouDF19UK4LPn2/eWG7G3GTKriks9luMczFnvnvqO15q8ZLlzmdKW/0B4H//U7ObK1Qwb1yapt2n0CeogwcPIoTA3T3Xa2ZZpA8+UAsQ/v67uv5kDoZkA9/6fssHuz8gKj6K/7X5H96e3pa7/MXt26p469y5asLtSy+p4eMFbZlhTSsiCn2COnToEE2aNClUAyT+/luN2nv9dbWigznsvbiXsdvHEnA9gK61uzK/53xcKriYJ5isxMTAwoX3qj88/7wal6+rP2iaRbPQPpi8IaXk8OHDhap7784d8PKCWrVUYyC/XYq6xOANg+m8vDORcZH8/NzP7By20zKTU2KiGolXrx5MmgStW6sLd+vW6eSkaQVAoW5BXbhwgfDwcNq0aWPuUPKElGqV8JAQtYK4k1P+nTveEM8XB75gxt8zSEpO4sNOH/LOY+9QzLZY/gWRXQkJsHKlKkUUEgIdOsCPP1reJDFN0zJVqBPUwYMHAQpNgpo/H376SfVUdeyYf+fdem4r438fz4VbF3i60dN80f0LapexwEnPCQmwYoVKTBcvgoeH6trr1UtXf9C0AqhQJ6hDhw5RvHhxXFwssPsph/75R62O26+f6q3KD+fDzzN+x3i2nd9GI+dG7HhhB93rds+fk+dEQoIaKv7xx2roeOvWaoRez546MWlaAVboE5SHhwfW1tbmDuWRXL+uVnWoWVP9Dpv6Nzc6IZoZe2fwxYEvcLBxYM4TcxjbZqzlVYFISIBly9SQ8dBQaNNGrdbYo4dOTJpWCBTaBBUXF4e/vz9vvfWWuUN5JAYDDBoEkZFqSHnp0qY7V7JMZu2JtUzaOYkrd64w3HU4s7rNopJTJdOdNDdiY1WmnjVLJaa2bWHxYujeXScmTStECm2COnbsGImJiQX++tN776kBEStWqILaprI7eDcT/5yI31U/y60CER6uVmJcsADCwqBdO1iyRC0frBOTphU6hXaY+aFDh4CCPUBiwQJVFm7kSFU1Ire893hn+NrpsNM8+cOTdF3ZlRt3b7Dy6ZUcHnHYspJTcLCq9lCjhppY27q1ytr//KNbTZpWiBXqBFW9enWqVKli7lByZc0a9Zvcr59KVI9ims+0h567eucqr/36Gs2+bsa+0H182u1Tzo45yzDXYZZToujoURg8WM1j+uYbNcH25Em1LlPnzjoxaVohV2i7+A4dOlRgW0/bt6vJuJ07q+k7Ntn8f8l7jzfent6Z7hOdEM2c/XOYs38OCUkJjG09limdpuBczPmR484TUsIff6im419/QcmSqoDruHFQtaq5o9M0LR8VygR148YNgoODeeONN8wdSo7t3w/PPquuN23ZAg4O2X/vNJ9pqQnKe4/3fS0nMU21NorbFudu4l2ea/IcMx+fSd2ydfMy/Ny7fVtdaFu0CM6eVesxzZ4Nr72mq4trWhFlIX05eaugXn86fhz69IFq1VQrKrPygZldVwLw9vRGTpUkTEkAoHZpNbG2ReUWHHzlID8995NlJKcTJ2DUKNU6GjdODVNcuVJdd5o4UScnTSvCCmUL6tChQ1hbW9OyZUtzh5Jtf/4JAwZAiRLqflYrP6S0ljJqKU3pOIXaZWozY+8MAMoVK8f8XvPpU78PwtzXbqKiVD28776Dw4fB3h6GDFGLBLZqZd7YNE2zGIU2Qbm6ulKsmAXWiUvH99+rkXqNG8Nvv0H16tl/r7fnvetOYpogenI03x/7nrkH5xISGYJ7FXfaVmvLmmfWmDcxJSfDvn1q/tK6darCuIuLWvpi2DDLWA5Y0zSLUugSVHJyMocPH2bo0KHmDiVLyclqXadPPlHFD3766eFuPc/lnuzx2gNkfF1pauepjHIfBUD1udW5FXeLDtU7sKDXAvO3mAIC1JDEtWvhv/9UhdshQ1TV29at9Ug8TdMyVOgSlK+vL1FRURZ//Sk6Wl3/X7tW/VYvWgS2tg/v53PRJ8tjbT6zmVn7ZgHgWcuTie0nmm8ek5TqutKmTbB+vRoWbm2tMvCsWWrcfPHi5olN07QCRUgpzR1Djri7u0tfX990XzMYDLRv356QkBDOnDlD2bJl8zm67NmxQy02ePGiaj29+27GDQkxTSCn3v//0Z34O5ScVRK3Sm74X/PHyc6JYc2HMb7teBqUa5APn+ABBgMcOgSbN6vEdOGC+kAdOqg6Tc8/D+XL539cmqaZlRDCT0qZ6+XMC1ULau7cuRw5coR169ZZZHIKD4c334RV/3nTyNGbffvUb/iDPJd73tdySunKK1+sPGExYanP+1/zB2CMxxhmdptp0tgfcvmyyrS//65GdURGqibg44+rcuv9+kHFivkbk6ZphUqhaUGdO3cOV1dXevXqxYYNG8w/Uo17E2cTE1VX3oQJcOsWGKYIYt+Rmc5xEtME8VPisZ9hz+utXmf96fWEx4ZTxqEMzzZ+lu+OfUfyh8n59zkvXYK9e9Xm4wNnzqjnq1RRy1r07Klq4pmymq2maQWKbkGhBka88sorODg4sHDhQpP8aGenSsODpvlMw/GQN199pX7fPTxg505w3ZTxBNxLUZfYfn47AOVmq5FtqwJW0bdhX4Y0HUKPej2ws7bju2PfmS453bmjygz5+sKRI2ooeHCweq1kSbVa4ssvq6TUtKke6KBpmkmYNEEJIXoC8wBr4Dsp5awHXrcHVgKtgHBgoJQyJKfnWbRoEfv27WPZsmVUrlw53X3SSzAZJZ30nk9bpSEzBoP6PV+9Gqiori89/jh0/MCbtVen4bpJ7ZfSbTey1UjcKrlx4NIBfj37KxFxEanHik6IBmBcm3F88vgn951nauepWcaSpbg4db3o5Ek4deredu6cGuwAahEqd3c1ibZTJ1XiooCvr6VpWsFgsi4+IYQ1cA54ArgEHAEGSylPp9lnNNBcSjlSCDEI6C+lHJjZcR/s4gsJCaFp06Y89thjbN++PcNWRXqDDdJ7Lqf7GgwQGKjKxv35J2yP8yah3cPFWce2Hkuver04FXaKiX9OpE3VNpwKO5WahCoUr0C7au14rMZj9KrXi6ZfN033fNkmJUREwJUrart8WW1BQfe2y5fvJSIrK1WU1cUF3NxUc69Vq6xnDGuapmXAkrv4WgOBUsogACHEj0A/4HSaffoB3sb764GvhBBCZpI1o6OiWblkBSEhwVy8GIqf7xFKO5bg3Tcn8a//KSSAlEgkMjkZabwP4H/YVz2WEimTATiy/wBSqv2QkkSDOs+f2/aw/L9l/HBtZeq5U1o9Le8MpuL5Z7h6M5brt+IxWMWQ7BBFiYqRtGwQRUn757AuHcH2639R3r4sYfERLDi8gAWH75UlL54o8KrRj9blmtG+nBt1ildXyTUpCf6LUzvt3g3x8Wrl2Pj4e1tcnKpdFxl5b7t16/7HkZEqez6oShWoW1c16+rUuZeUGjbMWeE/TdM0EzNlC2oA0FNK+arx8TCgjZRyTJp9Thr3uWR8fMG4z80Mj1tFSF43Sch5wt4AZWPBOQZOVITXfKFalNrqR0DjMFjQBrz3ZH4cb8+s98HRUQ1KSNnKlLn/cYUKKiFVrapuK1XSSUjTtHxjyS2oPCOEeA14DaBUyWI8c64jjvYO2NjYIBAgBAK4978CmXJPqMdzKq9hwpUXsBL39vy08kreueaFkCnvElhZC2ZU+B7viBFYW4GNlcDeAd6yW8x39iOxtxM42dtQzMoOR2GHo5UdxYQdpa2LU8aqOI7W9qlxe9/ahPfAZ+4NIjDeegO8dv9zqbdWVmBvj7e9PUy1V3Xq0ttKlVK3mqZphZQpW1DtAG8pZQ/j48kAUsqZafbZYdzngBDCBrgGlM+siy+zibqZxvOI16ByM4pP0zStKHvUFpQpl9s4AtQXQtQWQtgBg4AtD+yzBRhuvD8A2JVZcnoU6Y16y2gkXHrP6+SkaZqWv0w6UVcI0Rv4EjXMfKmU8mMhxHTAV0q5RQjhAKwCWgARwKCUQRUZyW0LStM0TctfFn0NSkq5Ddj2wHMfprkfBzxnyhg0TdO0gqlQrqiraZqmFXw6QWmapmkWSScoTdM0zSIVuGrmQogw4KK540iHM5DhBGMLVlDjBh27ORTUuEHHbg4NpZQlcvvmAjFRNy0ppUWufCeE8H2U0SrmUlDjBh27ORTUuEHHbg5CiEcacq27+DRN0zSLpBOUpmmaZpF0gso7i80dQC4V1LhBx24OBTVu0LGbwyPFXeAGSWiapmlFg25BaZqmaRZJJyhN0zTNIukElQtCiBAhxAkhhH/KMEohRFkhxJ9CiPPG2zLmjhNACLFUCHHDuDhkynPpxiqU+UKIQCFEgBCipfkizzB2byHEZeN3728sSJzy2mRj7GeFED3MEzUIIaoLIXYLIU4LIU4JIcYZn7fo7z2TuAvCd+4ghDgshDhujH2a8fnaQohDxhjXGVdWQAhhb3wcaHy9lgXGvlwIEZzme3czPm8R/17SxG8thDgmhNhqfJx33/m9JdD1lt0NCAGcH3huNvCu8f67wKfmjtMYSyegJXAyq1iB3sB2QABtgUMWGLs3MCGdfZsAxwF7oDZwAbA2U9yVgZbG+yWAc8b4LPp7zyTugvCdC8DJeN8WOGT8Ln9CrZIA8A0wynh/NPCN8f4gYJ054s4i9uXAgHT2t4h/L2nieQv4AdhqfJxn37luQeWdfsAK4/0VwNPmC+UeKeVe1FImaWUUaz9gpVQOAqWFEJXzJdB0ZBB7RvoBP0op46WUwUAg0NpkwWVCSnlVSnnUeP8O8C9QFQv/3jOJOyOW9J1LKWW08aGtcZNAV2C98fkHv/OU/y/WA48LkbKsdf7KJPaMWMS/FwAhRDWgD/Cd8bEgD79znaByRwJ/CCH8hFqOHqCilPKq8f41oKJ5QsuWjGKtCvyXZr9LZP4DZS5jjF0bS9N0pVpk7MZujBaov4oLzPf+QNxQAL5zY1eTP3AD+BPVoouUUhqMu6SNLzV24+u3gXL5GnAaD8YupUz53j82fu9zhRD2xucs6Xv/EpgEJBsflyMPv3OdoHLnMSllS6AX8IYQolPaF6VqwxaI8fsFKVajr4G6gBtwFfjcrNFkQgjhBGwAxkspo9K+ZsnfezpxF4jvXEqZJKV0A6qhWnKNzBtR9j0YuxCiKTAZ9Rk8gLLAO+aL8GFCiCeBG1JKP1OdQyeoXJBSXjbe3gA2of5juJ7SzDbe3jBfhFnKKNbLQPU0+1UzPmcxpJTXjf8xJwNLuNelZFGxCyFsUT/ya6SUG41PW/z3nl7cBeU7TyGljAR2A+1Q3V8pNUfTxpcau/H1UkB4/kb6sDSx9zR2uUopZTywDMv73jsAfYUQIcCPqK69eeThd64TVA4JIYoLIUqk3Ae6AyeBLcBw427DgV/ME2G2ZBTrFuBF4yihtsDtNF1SFuGBvvb+qO8eVOyDjCOFagP1gcP5HR+k9sN/D/wrpfwizUsW/b1nFHcB+c7LCyFKG+87Ak+grqHtBgYYd3vwO0/5/2IAsMvYqs13GcR+Js0fMwJ1HSft9272fy9SyslSympSylqoQQ+7pJRDycvvPL9GehSWDaiDGrl0HDgFvG98vhzwF3Ae2AmUNXesxrjWorplElH9wa9kFCtqVNBCVN/9CcDdAmNfZYwtwPgPvnKa/d83xn4W6GXGuB9Ddd8FAP7Grbelf++ZxF0QvvPmwDFjjCeBD43P10ElzUDgZ8De+LyD8XGg8fU6Fhj7LuP3fhJYzb2Rfhbx7+WBz+DJvVF8efad61JHmqZpmkXSXXyapmmaRdIJStM0TbNIOkFpmqZpFkknKE3TNM0i6QSlaZqmWSSdoLQiRQhRLk116Gvi/irdduaOLy0hhKcQon0+nStECOFsvN/KWEW7RX6cW9MyYpP1LppWeEgpw1ElexBCeAPRUso55opHCGEj79Ute5AnEA3sz6PjZef9zVGFPAdKKY/l9jialhd0C0or8owtBh9j8d8daWbw7zEW6fQVQvwrhPAQQmwUaj2nGcZ9agkhzggh1hj3WS+EKJaN434p1Fpi44QQTwm1Ps4xIcROIURFY7HWkcCbxtZdR6HWBxqQJu5o462nEOJvIcQW4LSx8OhnQogjxkKjr2fzq2gMbAaGSSnNUhFC09LSCUor6gSwALXuTitgKfBxmtcTpJTuqHVtfgHeAJoCXkKIlErMDYFFUsrGQBQw2ljTLrPj2kkp3aWUnwP7gLZSyhaommaTpJQhxnPOlVK6SSn/zuJztATGSSkboCpu3JZSeqAKjY4wliLKyi/AGCnlvmzsq2kmp7v4tKLOHpVw/lQlz7BGlVdKscV4ewI4JY01z4QQQajCl5HAf1LKf4z7rQb+B/yexXHXpblfDVhnbGHZAcG5+ByHpVqTCVR9yOZpWlulUHXysjruTuBVIcQOKWVSLmLQtDylE5RW1AlU4mmXwevxxtvkNPdTHqf89/NgvTCZjePeTXN/AfCFlHKLEMITtYJtegwYez2EEFaoZJbe8QQwVkq5I4PjZGQMqtW2CMhut6CmmYzu4tOKunigvBCiHajlJoQQLjk8Ro2U9wNDUF12Z3Nw3FLcW5JgeJrn76CWXk8RArQy3u+LWnk1PTuAUcZuRoQQDYyV9xFCnMnkcyQb428khJieyX6ali90gtKKumRU6f9PhRDHURW8czq0+yxq4cp/gTLA11LKhBwc1xv4WQjhB9xM8/yvQP+UQRKotZg6G4/XjvtbTWl9B5wGjgohTgLfAjbGYeSZLrEtpYxDJb++Qog3Mv/YmmZaupq5pj0C42i7rVLKpuaOJStCrYBaR0o539yxaFp26GtQmlZESCm3mjsGTcsJ3YLSNE3TLJK+BqVpmqZZJJ2gNE3TNIukE5SmaZpmkXSC0jRN0yySTlCapmmaRfo/Uve9vrRvY+UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "HeatCapacity = constant_volume_heat_capacity\n",
    "\n",
    "R = const.N_A * const.KB_J_PER_K\n",
    "plt.plot(temps,\n",
    "         HeatCapacity(eq, temps) / R, '-k',\n",
    "         label='Cv_rot Equilibrium')\n",
    "plt.plot(temps, HeatCapacity(para, temps) / R, '-b',\n",
    "         label='Cv_rot Para')\n",
    "plt.plot(temps, HeatCapacity(ortho, temps) / R, '-r',\n",
    "         label='Cv_rot Ortho')\n",
    "plt.plot(temps, 0.25 * HeatCapacity(para, temps) / R\n",
    "         + 0.75 * HeatCapacity(ortho, temps) / R, '-g',\n",
    "         label='Cv_rot 1:3 para:ortho')\n",
    "plt.plot(df_brink_T, df_brink_Cv, '+g')\n",
    "plt.plot(df_eucken_T, df_eucken_Cv, '+g')\n",
    "plt.plot(df_gia_T, df_gia_Cv, '+g')\n",
    "plt.plot(df_parting_T, df_parting_Cv, '+g')\n",
    "plt.plot(df_ce_T, df_ce_Cv, '+g', label = 'experimental data')\n",
    "plt.legend(loc='upper right', framealpha=100, frameon=False)\n",
    "plt.xlim(10, 400)\n",
    "plt.ylim(-0.1, 2.8)\n",
    "plt.xlabel('Temperature, K')\n",
    "plt.ylabel('Cv (rotational)/R')\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.23.0</td></tr><tr><td>Terra</td><td>0.16.0</td></tr><tr><td>Aer</td><td>0.7.0</td></tr><tr><td>Ignis</td><td>0.5.0</td></tr><tr><td>Aqua</td><td>0.8.0</td></tr><tr><td>IBM Q Provider</td><td>0.11.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.6.12 |Anaconda, Inc.| (default, Sep  8 2020, 17:50:39) \n",
       "[GCC Clang 10.0.0 ]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>6</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Wed Oct 21 18:45:59 2020 CEST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
